1. Experimental metadata

library(nlme)
library(aomisc)
## Warning: package 'carData' was built under R version 4.0.5
library(emmeans)
library(ggpubr)
library(tidyverse)
## Warning: package 'tidyr' was built under R version 4.0.5
## Warning: package 'readr' was built under R version 4.0.5
load(file = "1_data/psii_metadata.RData") 
psii = psii %>% mutate(Species = factor(Species, levels = c("A. cf humilis", "P. meandrina", "P. verrucosa")))
head(psii)
##          Region         Reef   Site       Date       Species     Depth
## 1 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis 10.967742
## 2 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis 10.967742
## 3 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis  8.387097
## 4 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis  8.387097
## 5 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis  8.387097
## 6 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis  8.387097
##   Bleaching catBleaching Treatment Tank Rack Clip      Vial max.yield
## 1         2         Poor        T0  2-2   24  530 A hum 275     0.471
## 2         2         Poor        T3  3-4   41  233 A hum 275     0.470
## 3         3         Fair        T0  2-2   24 1002 A hum 276     0.604
## 4         3         Fair        T3  3-4   41  265 A hum 276     0.624
## 5         3         Fair        T6  4-3   45  304 A hum 276     0.429
## 6         3         Fair        T9  1-4    5  992 A hum 276     0.235
##   mean.yield median.yield    sd.yield TargetTemp Temperature sd.Temperature
## 1  0.4346667        0.430 0.034239354      29.99       30.57           0.05
## 2  0.4700000        0.470 0.000000000      31.81       32.06           0.24
## 3  0.5773333        0.587 0.032593455      29.99       30.57           0.05
## 4  0.6203333        0.620 0.003511885      31.81       32.06           0.24
## 5  0.4290000        0.429          NA      34.81       34.97           0.36
## 6  0.1950000        0.219 0.056000000      37.81       36.33           0.90
##     MMM MMM_dev DHW duration int_mean  int_max  int_cum
## 1 28.96    1.61  10       32 1.386897 2.007517 12.48207
## 2 28.96    3.10  10       32 1.386897 2.007517 12.48207
## 3 28.96    1.61  10       32 1.386897 2.007517 12.48207
## 4 28.96    3.10  10       32 1.386897 2.007517 12.48207
## 5 28.96    6.01  10       32 1.386897 2.007517 12.48207
## 6 28.96    7.87  10       32 1.386897 2.007517 12.48207

We’ve combined data from the sampling of coral nubbins, experimental data and results, with long-term environmental data into a single data frame for analysis. The data include the following information:

Geographic metadata Region: Bioregions reflective of the biogeoraphy of coral and fish diversity in the Coral Sea Marine Park Reef: Distinct coral reef groups Site: Sampling site

Sampling metadata Date: Date of collection Species: Species name. Pocillopora confirmed using RFLP Depth: Collection depth Bleaching: Bleaching category 1 (bleached) to 6 (healthy) catBleaching: Bleaching category grouping (Poor: 1&2, Fair: 3&4, Good:5&6)

Experimental metadata Treatment: Temperature treatment Tank: Tank within treatment jackets, 3 tanks per jacket Rack: Rack within Tank Clip: Clip on rack Vial: Unique identification number for each coral genet

Experimental measurements max.yield: maximum measurement of PSII activity from 3 measures mean.yield: mean measurement of PSII activity from 3 measures median.yield: median measurement of PSII activity from 3 measures sd.yield: standard deviation of PSII activity from 3 measures

Treatment data TargetTemp: Target treatment temperature Temperature: Average measured treatment temperature sd.Temperature: Standard deviation of measured treatment temperature

NOAA Coral Reef Watch metadata MMM: Local Maximum Monthly Mean sea surface temperature at the time of sampling MMM_dev: Deviation from MM at the time of sampling DHW: Degree Heating Weeks at the time of sampling

RmarineHeatWaves metadata (based on NOAA CRW SST data) duration: Heatwave duration (days) int_mean: Mean intensity int_max: Max intensity int_cum: Cumulative intensity

2. Dose Response curves

In order to run a log logistic regression on the experimental data and identify differences between species we need to make sure all species have a common set of reefs and enough samples from each reef and bleaching categories to explore the effects of any interactions. So we’ll create two new data frames, one for model exploration/selection (psii.1) and another for the comparison between species (psii.2). These exclude some reefs where sample sizes were low.

#How many samples from each Species/Reef. This will include all fragments (sum of all fragments, NOT individual genotypes)
psii %>% group_by(Species, Reef) %>% 
  dplyr::summarise(n = dplyr::n()) %>% 
  spread(key = Species, value = n)
## `summarise()` has grouped output by 'Species'. You can override using the
## `.groups` argument.
## # A tibble: 9 × 4
##   Reef         `A. cf humilis` `P. meandrina` `P. verrucosa`
##   <fct>                  <int>          <int>          <int>
## 1 Bougainville              73             36             59
## 2 Chilcott                  70             39             11
## 3 Flinders                 112             20             97
## 4 Frederick                 62             48              6
## 5 Herald                    62             16             36
## 6 Lihou                    117             23             71
## 7 Moore                     67             63             37
## 8 Saumarez                  NA             85             NA
## 9 Wreck                     96             52             24
#Creating a subset of the data for model selection
psii.1 = psii %>% filter(! Reef %in% c("Saumarez", "Frederick", "Chilcott", "Wreck")) %>%
  mutate(Reef = factor(Reef)) %>%
  mutate(Species = factor(Species)) %>% 
  droplevels() %>% 
  as.data.frame()

#Creating a subset of the data for model selection
psii.2 = psii %>% filter(! Reef %in% c("Saumarez", "Frederick")) %>%
  mutate(Reef = factor(Reef)) %>%
  mutate(Species = factor(Species)) %>% 
  droplevels() %>% 
  as.data.frame()

#Table 1. Number of individual corals collected from each reef AND used in the model.

#count of how many individuals were collected per reef 
table1 = psii %>% group_by(Species, Reef, Vial) %>% 
dplyr::summarise() %>% 
  group_by(Species, Reef) %>% 
  dplyr::summarise(ln = n()) %>% 
  spread(Species, ln) 
## `summarise()` has grouped output by 'Species', 'Reef'. You can override using
## the `.groups` argument.
## `summarise()` has grouped output by 'Species'. You can override using the
## `.groups` argument.
table1
## # A tibble: 9 × 4
##   Reef         `A. cf humilis` `P. meandrina` `P. verrucosa`
##   <fct>                  <int>          <int>          <int>
## 1 Bougainville              23             10             17
## 2 Chilcott                  18             10              3
## 3 Flinders                  30              5             26
## 4 Frederick                 17             13              2
## 5 Herald                    18              5             11
## 6 Lihou                     32              6             18
## 7 Moore                     20             17             10
## 8 Saumarez                  NA             22             NA
## 9 Wreck                     24             13              6
psii %>% group_by(Species, Reef, Treatment) %>% 
  summarise(mean = round(mean(median.yield),2)) %>% 
  spread(key = Treatment, value = mean)
## `summarise()` has grouped output by 'Species', 'Reef'. You can override using
## the `.groups` argument.
## # A tibble: 25 × 6
## # Groups:   Species, Reef [25]
##    Species       Reef            T0    T3    T6    T9
##    <fct>         <fct>        <dbl> <dbl> <dbl> <dbl>
##  1 A. cf humilis Bougainville  0.51  0.55  0.45  0.1 
##  2 A. cf humilis Chilcott      0.62  0.63  0.47  0.07
##  3 A. cf humilis Flinders      0.61  0.61  0.44  0.08
##  4 A. cf humilis Frederick     0.62  0.62  0.51  0.08
##  5 A. cf humilis Herald        0.6   0.59  0.35  0.04
##  6 A. cf humilis Lihou         0.61  0.61  0.49  0.03
##  7 A. cf humilis Moore         0.55  0.57  0.41  0.07
##  8 A. cf humilis Wreck         0.65  0.65  0.61  0.18
##  9 P. meandrina  Bougainville  0.61  0.61  0.58  0.1 
## 10 P. meandrina  Chilcott      0.64  0.65  0.56  0.05
## # … with 15 more rows
## # ℹ Use `print(n = ...)` to see more rows
psii %>% group_by(Species, Reef, Treatment) %>% 
  summarise(mean = round(sd(median.yield),2)) %>% 
  spread(key = Treatment, value = mean)
## `summarise()` has grouped output by 'Species', 'Reef'. You can override using
## the `.groups` argument.
## # A tibble: 25 × 6
## # Groups:   Species, Reef [25]
##    Species       Reef            T0    T3    T6    T9
##    <fct>         <fct>        <dbl> <dbl> <dbl> <dbl>
##  1 A. cf humilis Bougainville  0.08  0.06  0.07  0.07
##  2 A. cf humilis Chilcott      0.03  0.02  0.08  0.05
##  3 A. cf humilis Flinders      0.04  0.03  0.13  0.1 
##  4 A. cf humilis Frederick     0.02  0.03  0.18  0.08
##  5 A. cf humilis Herald        0.03  0.04  0.09  0.04
##  6 A. cf humilis Lihou         0.03  0.03  0.13  0.02
##  7 A. cf humilis Moore         0.06  0.05  0.13  0.08
##  8 A. cf humilis Wreck         0.02  0.02  0.05  0.11
##  9 P. meandrina  Bougainville  0.08  0.07  0.07  0.08
## 10 P. meandrina  Chilcott      0.02  0.04  0.07  0.04
## # … with 15 more rows
## # ℹ Use `print(n = ...)` to see more rows

Results 3.2 Differences in Fv/Fm observed in different treatments

#reporting differences among treatments in results section 3.2
psii.2 %>% 
  dplyr::group_by(Treatment) %>% 
  summarise(med.treat = mean(median.yield), sd.treat = sd(median.yield), median.treat = median(median.yield))
## # A tibble: 4 × 4
##   Treatment med.treat sd.treat median.treat
##   <fct>         <dbl>    <dbl>        <dbl>
## 1 T0           0.610    0.0567        0.624
## 2 T3           0.619    0.0486        0.628
## 3 T6           0.546    0.114         0.588
## 4 T9           0.0855   0.0899        0.05

Assuming a fixed lower asymptote representing complete loss of photosynthetic yield, a three parameter log-logistic model can be assumed, estimating the upper asymptote, the location of the inflection point (PSII-50), and steepness of the curve whereby:

b is the slope around the inflection point d is the upper asymptote, e is the inflection point (PSII-50)

2.1 Model selection : SPECIES

What parameters to estimate?

Our primary interest is in the inflection point (e or ED50) but other parameter estimates may be important and may vary among species. We anticipate most of the variation in a species response to temperature treatment at this scale will occur at the Reef level so we included Reef as a random factor to explore the importance of parameter estimates

#Naive model to identify starting points
# Model Parameters:  b = shape parameter (i.e steepness); d=upper limit; e=ED50
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species,
            pmodels = c( ~ 1,  ~ 1,  ~ Species), bcVal = 0.5)

species.1 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef,
                    fixed = list(b ~ 1, d ~ 1, e ~ Species),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species,
            pmodels = c( ~ 1,  ~ Species,  ~ Species), bcVal = 0.5)

species.2 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef,
                    fixed = list(b ~ 1, d ~ Species, e ~ Species),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species,
            pmodels = c( ~ Species,  ~ 1,  ~ Species), bcVal = 0.5)

species.3 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef,
                    fixed = list(b ~ Species, d ~ 1, e ~ Species),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species,
            pmodels = c( ~ Species,  ~ Species,  ~ Species), bcVal = 0.5)

species.4 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef,
                    fixed = list(b ~ Species, d ~ Species, e ~ Species),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

AIC(species.1,species.2,species.3, species.4)
##           df       AIC
## species.1  7 -2131.489
## species.2  9 -2179.997
## species.3  9 -2171.099
## species.4 11 -2203.347

The more complex model that estimates all 3 parameters (species.4) is the best model. However, it may be necessary to reduce the complexity of the model when accounting for other random factors in the model, in which case estimating d and e provide the best result (species.2). For example, when accounting for Depth and Tank effects, the models are too complex to also estimate b.

Accounting for Depth effects

modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species,
            pmodels = c( ~ 1,  ~ Species,  ~ Species), bcVal = 0.5)

species.5 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef/Depth,
                    fixed = list(b ~ 1, d ~ Species, e ~ Species),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

AIC(species.2,species.5) 
##           df       AIC
## species.2  9 -2179.997
## species.5 10 -2177.997
anova(species.2, species.5)
##           Model df       AIC       BIC   logLik   Test      L.Ratio p-value
## species.2     1  9 -2179.997 -2136.886 1098.998                            
## species.5     2 10 -2177.997 -2130.096 1098.998 1 vs 2 8.668175e-06  0.9977

Depth does not need to be accounted for.

Accounting for Tank effects

modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species,
            pmodels = c( ~ 1,  ~ Species,  ~ Species), bcVal = 0.5)

species.6 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef/Tank,
                    fixed = list(b ~ 1, d ~ Species, e ~ Species),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(species.2,species.6)
##           df       AIC
## species.2  9 -2179.997
## species.6 10 -2177.081
anova(species.2, species.6)
##           Model df       AIC       BIC   logLik   Test   L.Ratio p-value
## species.2     1  9 -2179.997 -2136.886 1098.998                         
## species.6     2 10 -2177.081 -2129.180 1098.541 1 vs 2 0.9153679  0.3387

Tank does not need to be accounted for.

Accounting for bleaching

Introduce an interaction between Species and Bleaching categories. We explore the interaction only for d and e.

#effect of bleached samples on e
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species:catBleaching,
            pmodels = c( ~ 1,  ~ Species,  ~ Species*catBleaching), bcVal = 0.5)

species.7 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef,
                    fixed = list(b ~ 1, d ~ Species, e ~ Species*catBleaching),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

#effect of bleached samples on d
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species:catBleaching,
            pmodels = c( ~ 1,  ~ Species*catBleaching,  ~ Species), bcVal = 0.5)

species.8 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef,
                    fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

#effect of bleached samples on d and e
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.1,
            curveid = Species:catBleaching,
            pmodels = c( ~ 1,  ~ Species*catBleaching,  ~ Species*catBleaching), bcVal = 0.5)

species.9 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = e ~ 1|Reef,
                    fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species*catBleaching),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

#check random factors on d and e
species.9.1 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
                    random = d + e ~ 1|Reef,
                    fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species*catBleaching),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.1, : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT
## message: function evaluation limit reached without convergence (9)
AIC(species.2, species.7, species.8, species.9, species.9.1)
##             df       AIC
## species.2    9 -2179.997
## species.7   15 -2181.960
## species.8   15 -2230.468
## species.9   21 -2244.262
## species.9.1 23 -2277.888
anova(species.2,species.9.1)
##             Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## species.2       1  9 -2179.997 -2136.886 1098.998                        
## species.9.1     2 23 -2277.888 -2167.716 1161.944 1 vs 2 125.8915  <.0001
anova(species.9.1) # interaction is significant
##                        numDF denDF   F-value p-value
## b                          1   866   875.517  <.0001
## d.(Intercept)              1   866   315.787  <.0001
## d.Species                  2   866    26.764  <.0001
## d.catBleaching             2   866   117.710  <.0001
## d.Species:catBleaching     4   866     9.038  <.0001
## e.(Intercept)              1   866 24674.744  <.0001
## e.Species                  2   866    89.599  <.0001
## e.catBleaching             2   866    10.094  <.0001
## e.Species:catBleaching     4   866     2.932    0.02

The interaction between bleaching categories and species is important and improves the model so we need to control for the level of bleaching of sampled nubbins. The effect on e and d is small but should not be ignored. If we look at the contrast below, the species differences change depending on the bleaching category.

#The effect of including the interaction between species and catBleaching on d and e is small 
emmeans(species.9.1,pairwise~catBleaching|Species, param="d", level = 0.95)$emmeans %>% plot(comparisons = TRUE)

emmeans(species.9.1,pairwise~catBleaching|Species, param="e", level = 0.95)$emmeans %>% plot(comparisons = TRUE)

#The contrast show how bleaching can affect the species comparisons
emmeans(species.9.1,pairwise~catBleaching|Species, param="d", level = 0.95)$contrast %>% plot()

emmeans(species.9.1,pairwise~Species|catBleaching, param="d", level = 0.95)$contrast %>% plot()

emmeans(species.9.1,pairwise~catBleaching|Species, param="e", level = 0.95)$contrast %>% plot()

emmeans(species.9.1,pairwise~Species|catBleaching, param="e", level = 0.95)$contrast %>% plot()

Note: The blue bars are confidence intervals for the EMMs, and the red arrows are for the comparisons among them. If an arrow from one mean overlaps an arrow from another group, the difference is not “significant,” based on the adjust setting (which defaults to “tukey”) and the value of alpha (which defaults to 0.05). (https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html).

2.2 Model predictions: Species*Bleaching

#including larger subset of reefs (psii.2)
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.2,
            curveid = Species:catBleaching,
            pmodels = c( ~ 1,  ~ Species*catBleaching,  ~ Species*catBleaching), bcVal = 0.5)

species.9.2 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
                    random = d + e ~ 1|Reef,
                    fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species*catBleaching),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT
## message: function evaluation limit reached without convergence (9)
#Chilcott and Wreck have a strong effect on parameter estimates of e and d for each species particularly for A. cf humilis
emmeans(species.9.1, ~Species, param = 'e')
## NOTE: Results may be misleading due to involvement in interactions
##  Species       emmean     SE  df lower.CL upper.CL
##  A. cf humilis   6.76 0.0607 866     6.64     6.87
##  P. meandrina    7.33 0.0814 866     7.17     7.49
##  P. verrucosa    7.73 0.0663 866     7.60     7.86
## 
## Results are averaged over the levels of: catBleaching 
## Confidence level used: 0.95
emmeans(species.9.2, ~Species, param = 'e')
## NOTE: Results may be misleading due to involvement in interactions
##  Species       emmean    SE   df lower.CL upper.CL
##  A. cf humilis   7.10 0.140 1156     6.82     7.37
##  P. meandrina    7.36 0.143 1156     7.08     7.64
##  P. verrucosa    7.73 0.143 1156     7.45     8.01
## 
## Results are averaged over the levels of: catBleaching 
## Confidence level used: 0.95
emmeans(species.9.1, ~Species, param = 'd')
## NOTE: Results may be misleading due to involvement in interactions
##  Species       emmean      SE  df lower.CL upper.CL
##  A. cf humilis  0.578 0.00814 866    0.562    0.594
##  P. meandrina   0.624 0.00970 866    0.605    0.644
##  P. verrucosa   0.618 0.00845 866    0.602    0.635
## 
## Results are averaged over the levels of: catBleaching 
## Confidence level used: 0.95
emmeans(species.9.2, ~Species, param = 'd')
## NOTE: Results may be misleading due to involvement in interactions
##  Species       emmean      SE   df lower.CL upper.CL
##  A. cf humilis  0.585 0.00614 1156    0.573    0.597
##  P. meandrina   0.632 0.00747 1156    0.618    0.647
##  P. verrucosa   0.628 0.00692 1156    0.615    0.642
## 
## Results are averaged over the levels of: catBleaching 
## Confidence level used: 0.95
#what about the slope b
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.2,
            curveid = Species:catBleaching,
            pmodels = c( ~ Species*catBleaching,  ~ Species*catBleaching,  ~ Species*catBleaching), bcVal = 0.5)

species.9.3 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
                    random = b + d + e ~ 1|Reef,
                    fixed = list(b ~ Species*catBleaching, d ~ Species*catBleaching, e ~ Species*catBleaching),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT
## message: function evaluation limit reached without convergence (9)
## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Singular precision matrix in level -1, block 1

## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Singular precision matrix in level -1, block 1

## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Singular precision matrix in level -1, block 1
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.2,
            curveid = Species:catBleaching,
            pmodels = c( ~ Species,  ~ Species*catBleaching,  ~ Species*catBleaching), bcVal = 0.5)

#Does not run
#species.9.4 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
#                    random = b + d + e ~ 1|Reef,
#                    fixed = list(b ~ Species, d ~ Species*catBleaching, e ~ Species*catBleaching),
#                    start = coef(modNaive.Species), 
#                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.2,
            curveid = Species:catBleaching,
            pmodels = c( ~ catBleaching,  ~ Species*catBleaching,  ~ Species*catBleaching), bcVal = 0.5)

species.9.5 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
                    random = b + d + e ~ 1|Reef,
                    fixed = list(b ~ catBleaching, d ~ Species*catBleaching, e ~ Species*catBleaching),
                    start = coef(modNaive.Species), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

AIC(species.9.2, species.9.3, species.9.5)
##             df       AIC
## species.9.2 23 -2982.373
## species.9.3 34 -3109.717
## species.9.5 28 -3034.841
anova(species.9.2, species.9.3)
##             Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## species.9.2     1 23 -2982.373 -2865.669 1514.187                        
## species.9.3     2 34 -3109.718 -2937.198 1588.859 1 vs 2 149.3442  <.0001

Species also have different slopes and so it is important that we consider this in our model. Significant improvements to model fit with species.9.3

Estimated marginal mean & pairwise contrasts

parameter e

my.model1 = species.9.3
load("3_outputs/nlme_All species.RData")
my.data1 = psii.2

emmeans(my.model1, ~Species, param = 'e', level = 0.95) %>% plot(comparisons = TRUE)
## NOTE: Results may be misleading due to involvement in interactions

emmeans(my.model1, pairwise~Species, param = 'e', level = 0.95)$contrast %>% plot()
## NOTE: Results may be misleading due to involvement in interactions

The differences in the inflection point (e, ED50) between species are significant.

emmeans(my.model1, pairwise~Species, param = 'e', level = 0.95)$contrast    #ED50 contrasts 
## NOTE: Results may be misleading due to involvement in interactions
##  contrast                     estimate     SE   df t.ratio p.value
##  A. cf humilis - P. meandrina   -0.376 0.0765 1148 -4.914  <.0001 
##  A. cf humilis - P. verrucosa   -0.696 0.0767 1148 -9.067  <.0001 
##  P. meandrina - P. verrucosa    -0.320 0.0857 1148 -3.733  0.0006 
## 
## Results are averaged over the levels of: catBleaching 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans(my.model1, ~Species, param = 'e', level = 0.95) 
## NOTE: Results may be misleading due to involvement in interactions
##  Species       emmean    SE   df lower.CL upper.CL
##  A. cf humilis   7.05 0.154 1148     6.75     7.35
##  P. meandrina    7.42 0.159 1148     7.11     7.74
##  P. verrucosa    7.74 0.160 1148     7.43     8.06
## 
## Results are averaged over the levels of: catBleaching 
## Confidence level used: 0.95

parameter d

emmeans(my.model1, ~Species, param = 'd', level = 0.95) %>% plot(comparisons = TRUE)
## NOTE: Results may be misleading due to involvement in interactions

emmeans(my.model1, pairwise~Species, param = 'd', level = 0.95)$contrast %>% plot()
## NOTE: Results may be misleading due to involvement in interactions

A. cf humilis has lower upper asymptote than either P. meandrina or P. verrucosa. There is no observable difference between Pocilloporids.

emmeans(my.model1, pairwise~Species, param = 'd', level = 0.95)$contrast    #ED50 contrasts 
## NOTE: Results may be misleading due to involvement in interactions
##  contrast                     estimate      SE   df t.ratio p.value
##  A. cf humilis - P. meandrina -0.03078 0.00696 1148 -4.421  <.0001 
##  A. cf humilis - P. verrucosa -0.03531 0.00657 1148 -5.378  <.0001 
##  P. meandrina - P. verrucosa  -0.00453 0.00781 1148 -0.581  0.8306 
## 
## Results are averaged over the levels of: catBleaching 
## P value adjustment: tukey method for comparing a family of 3 estimates

parameter b

emmeans(my.model1, ~Species, param = 'b', level = 0.95) %>% plot(comparisons = TRUE)
## NOTE: Results may be misleading due to involvement in interactions

emmeans(my.model1, pairwise~Species, param = 'b', level = 0.95)$contrast %>% plot()
## NOTE: Results may be misleading due to involvement in interactions

A. cf humilis has slacker slope than either P. meandrina or P. verrucosa. The condition of A. cf humilis may have dropped at lower temps. There is no observable difference between Pocilloporids.

emmeans(my.model1, pairwise~Species, param = 'b', level = 0.95)$contrast    #ED50 contrasts 
## NOTE: Results may be misleading due to involvement in interactions
##  contrast                     estimate    SE   df t.ratio p.value
##  A. cf humilis - P. meandrina    0.853 0.156 1148  5.455  <.0001 
##  A. cf humilis - P. verrucosa    0.663 0.178 1148  3.722  0.0006 
##  P. meandrina - P. verrucosa    -0.190 0.218 1148 -0.875  0.6565 
## 
## Results are averaged over the levels of: catBleaching 
## P value adjustment: tukey method for comparing a family of 3 estimates

catBleaching interaction

#Including bleaching
emmeans(my.model1, ~catBleaching|Species, param = 'e', level = 0.95) %>% plot()

emmeans(my.model1, ~catBleaching|Species, param = 'd', level = 0.95) %>% plot()

emmeans(my.model1, ~catBleaching|Species, param = 'b', level = 0.95) %>% plot()

catBleaching has its strongest effect on d within species, whereby colonies that started in a poorer condition have lower values of median PSII yield.

SOM Tables S4 and S4: Species Pairwise Contrasts and Emmeans

emmeans(my.model1, pairwise~Species, param = 'e', level = 0.95)$contrast    #ED50 contrasts 
## NOTE: Results may be misleading due to involvement in interactions
##  contrast                     estimate     SE   df t.ratio p.value
##  A. cf humilis - P. meandrina   -0.376 0.0765 1148 -4.914  <.0001 
##  A. cf humilis - P. verrucosa   -0.696 0.0767 1148 -9.067  <.0001 
##  P. meandrina - P. verrucosa    -0.320 0.0857 1148 -3.733  0.0006 
## 
## Results are averaged over the levels of: catBleaching 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans(my.model1, pairwise~Species, param = 'd', level = 0.95)$contrast    #ED50 contrasts 
## NOTE: Results may be misleading due to involvement in interactions
##  contrast                     estimate      SE   df t.ratio p.value
##  A. cf humilis - P. meandrina -0.03078 0.00696 1148 -4.421  <.0001 
##  A. cf humilis - P. verrucosa -0.03531 0.00657 1148 -5.378  <.0001 
##  P. meandrina - P. verrucosa  -0.00453 0.00781 1148 -0.581  0.8306 
## 
## Results are averaged over the levels of: catBleaching 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans(my.model1, pairwise~Species, param = 'b', level = 0.95)$contrast    #ED50 contrasts 
## NOTE: Results may be misleading due to involvement in interactions
##  contrast                     estimate    SE   df t.ratio p.value
##  A. cf humilis - P. meandrina    0.853 0.156 1148  5.455  <.0001 
##  A. cf humilis - P. verrucosa    0.663 0.178 1148  3.722  0.0006 
##  P. meandrina - P. verrucosa    -0.190 0.218 1148 -0.875  0.6565 
## 
## Results are averaged over the levels of: catBleaching 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans(my.model1, ~Species, param = 'e', level = 0.95) 
## NOTE: Results may be misleading due to involvement in interactions
##  Species       emmean    SE   df lower.CL upper.CL
##  A. cf humilis   7.05 0.154 1148     6.75     7.35
##  P. meandrina    7.42 0.159 1148     7.11     7.74
##  P. verrucosa    7.74 0.160 1148     7.43     8.06
## 
## Results are averaged over the levels of: catBleaching 
## Confidence level used: 0.95

Model check

resid(my.model1) %>% plot()

VarCorr(my.model1)
## Reef = pdLogChol(list(b ~ 1,d ~ 1,e ~ 1)) 
##               Variance     StdDev       Corr         
## b.(Intercept) 2.979937e-20 1.726249e-10 b.(In) d.(In)
## d.(Intercept) 4.497247e-04 2.120671e-02 -0.662       
## e.(Intercept) 1.514108e-01 3.891154e-01 -0.935  0.374
## Residual      3.829744e-03 6.188493e-02
ranef(my.model1)
##              b.(Intercept) d.(Intercept) e.(Intercept)
## Bougainville  1.539353e-10  -0.036360850  -0.130969244
## Chilcott      3.784724e-11   0.022342256  -0.293119010
## Flinders     -7.389782e-11  -0.003692153   0.239259916
## Herald        1.733321e-10  -0.001971101  -0.472525825
## Lihou        -1.780246e-11   0.006769950  -0.006246521
## Moore         9.070098e-11  -0.015226812  -0.128861706
## Wreck        -3.641153e-10   0.028138711   0.792462390
#coef(my.model1)
#vcov(my.model1)

2.3 Figure: Species comparison

##Create predictions data frame
predframe.spp = with(my.data1, 
                     expand.grid(Species = levels(Species),
                                 catBleaching = levels(catBleaching),
                                 MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))

#Adding confidence intervals
prdc.spp = as.data.frame(nlraa::predict_nlme(my.model1, newdata = predframe.spp, 
                                           interval = "confidence", level = 0.95))

prdc.spp$Species = predframe.spp$Species
prdc.spp$MMM_dev = predframe.spp$MMM_dev
prdc.spp$catBleaching = predframe.spp$catBleaching

#Averaging predictions across catBleaching
prdc.average = prdc.spp %>% dplyr::select(-catBleaching) %>%
  group_by(Species, MMM_dev) %>% summarise_all(.funs=mean)

# Predict ED50 for each species
ED50spp = emmeans(my.model1, ~Species|catBleaching, param= "e") %>% as_tibble()

predframe.ED50 = ED50spp %>% dplyr::select(Species, catBleaching, emmean) %>% 
  dplyr::rename(MMM_dev = emmean)
  
prdc.ED50 = as.data.frame(nlraa::predict_nlme(my.model1, newdata = predframe.ED50, 
                                           interval = "confidence", level = 0.95))
prdc.ED50 = bind_cols(predframe.ED50, prdc.ED50)

prdc.ED50.average = prdc.ED50 %>% dplyr::select(-catBleaching) %>%
  group_by(Species) %>% summarise_all(.funs=mean)

Fig 2A: PSII comparison by species

library(wesanderson)
col = wes_palette("GrandBudapest1", 4, type = "continuous")

Species.plot = ggplot() +
  geom_segment(data = prdc.ED50.average, aes(x = MMM_dev, xend = MMM_dev, y = -Inf, 
                                             yend = Estimate , color = Species), 
                  linetype = 1, size = .3) +
  geom_point(data=psii.1, aes(x=MMM_dev, y=median.yield, col=Species), alpha=0.1, size = 1) +
  geom_line(data = prdc.average, aes(colour=Species, y=Estimate, x=MMM_dev), size = .3) +
  geom_ribbon(data=prdc.average, aes(x = MMM_dev, ymin = Q2.5, ymax = Q97.5, fill = Species), 
              alpha = 0.5, show.legend = FALSE) +
  geom_point(data = prdc.ED50.average, aes(x = MMM_dev, y = Estimate, fill = Species),
             size = 1.5, stroke = .2, shape = 21, col = "white") +
  scale_x_continuous(breaks = c(3,6,9)) +
  coord_cartesian(xlim = c(2,9.6), ylim = c(0,.75)) +
  scale_fill_manual(values=col) + 
  scale_colour_manual(values=col) +
  labs(x="Temp. above local MMM (°C)", 
       y="F[v]/F[m]") +
  theme_classic() +
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        legend.title = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        #legend.position = "none",
        legend.text = element_text(face = "italic", size = 7, family = "Helvetica"),
        legend.key.size = unit(2, units = "mm"),
        legend.margin = margin(0,0,0,0, unit = "mm"),
        axis.text.y = element_text(size = 8, family = "Helvetica"),
        axis.title = element_text(size = 9, family = "Helvetica"))  

Species.plot

ggsave("2_figures/SpeciesPlotonly.png", width = 140, height = 70, units = "mm")

Fig SX: PSII Species X catBleaching

SpeciesBl.plot = ggplot() +
  geom_segment(data = prdc.ED50, aes(x = MMM_dev, xend = MMM_dev, y = -Inf, yend = Estimate , color = Species), 
                  linetype = 1, size = .3) +
  geom_point(data=psii.1, aes(x=MMM_dev, y=median.yield, col=Species), alpha=0.1, size = 1) +
  geom_line(data = prdc.spp, aes(colour=Species, y=Estimate, x=MMM_dev), size = .3) +
  geom_ribbon(data=prdc.spp, aes(x = MMM_dev, ymin = Q2.5, ymax = Q97.5, fill = Species), 
              alpha = 0.5, show.legend = FALSE) +
  geom_point(data = prdc.ED50, aes(x = MMM_dev, y = Estimate, fill = Species), size = 1.5, stroke = .2, shape = 21, col = "white") +
  scale_x_continuous(breaks = c(3,6,9)) +
  coord_cartesian(xlim = c(2,9.6), ylim = c(0,.75)) +
  scale_fill_manual(values=col) + 
  scale_colour_manual(values=col) +
  labs(x="Temp. above local MMM (°C)", 
       y="Fv/Fm") +
  theme_classic() +
  facet_wrap(~catBleaching) +
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        legend.title = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        legend.position = "right",
        legend.text = element_text(face = "italic", size = 7, family = "Helvetica"),
        legend.key.size = unit(2, units = "mm"),
        legend.margin = margin(0,0,0,0, unit = "mm"),
        axis.text.y = element_text(size = 8, family = "Helvetica"),
        axis.title = element_text(size = 9, family = "Helvetica"))  

SpeciesBl.plot

ggsave("2_figures/FigS3_catBleaching_by_Species.png", width = 140, height = 70, units = "mm")

2.4 Model selection : REEF

Combined 3-way interaction

Since we don’t have the same set of reefs for all species, we can’t run the three-way interaction Species:Reef:catBleaching but we can still explore what this will look like for a subset of reefs and how it affects each model parameter. However, this significantly increases the complexity of the model so we need to simplify it.

If we look at the interaction between species and catBleaching. We saw that for species, the most complex model, with an interaction between species and catBleaching on all parameter estimates was important. However, the strength of the effect depends on the species and parameter. We can see from below: - catBleaching has a weak effect on b, though can be an important source of variation between species - catBleaching has a strong effect on d, particularly for A. humilis - catBleaching has a weak effect on e, with the only significant effect on P. verrucosa

emmeans(my.model1, ~catBleaching|Species, param = 'b', level = 0.95) %>% plot(comparisons = TRUE)

emmeans(my.model1, ~catBleaching|Species, param = 'd', level = 0.95) %>% plot(comparisons = TRUE)

emmeans(my.model1, ~catBleaching|Species, param = 'e', level = 0.95) %>% plot(comparisons = TRUE)

Therefore when developing a more complex model that include reefs. We must consider: - Species differences on b - Species and catBleahching differences on d - Species differences on e

However, running Species on b does not work, perhaps because it is confounded by another fixed factor. So we have to fix the slope.

#testing the 3-interaction on psii.2 (a common subset of reefs)
modNaive <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = psii.2,
            curveid = Species:Reef:catBleaching,
            pmodels = c( ~ 1,  ~ Species*catBleaching,  ~ Species*Reef), bcVal = 0.5)

reef.1 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
                    random = e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species*Reef),
                    start = coef(modNaive), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

emmeans(reef.1, ~Reef|Species, param = "e") %>% plot(comparisons = TRUE)
## Warning: Comparison discrepancy in group "P. verrucosa", Bougainville - Wreck:
##     Target overlap = 0.0101, overlap on graph = -0.0223

emmeans(reef.1, ~catBleaching|Species, param = "d") %>% plot(comparisons = TRUE)

anova(reef.1)
##                        numDF denDF   F-value p-value
## b                          1   829   8603.21  <.0001
## d.(Intercept)              1   829 151066.39  <.0001
## d.Species                  2   829     58.34  <.0001
## d.catBleaching             2   829    191.14  <.0001
## d.Species:catBleaching     4   829     30.84  <.0001
## e.(Intercept)              1   829 100928.49  <.0001
## e.Species                  2   829     13.86  <.0001
## e.Reef                     6   829     85.86  <.0001
## e.Species:Reef            12   829     11.11  <.0001
##Create predictions data frame
predframe.spp = with(psii.2, expand.grid(Species = levels(Species),
                                         Reef = levels(Reef),
                                         catBleaching = levels(catBleaching),
                                         MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))

#Adding confidence intervals
prdc.spp = as.data.frame(nlraa::predict_nlme(reef.1, newdata = predframe.spp, 
                                           interval = "confidence", level = 0.950))

prdc.spp = bind_cols(predframe.spp, prdc.spp)

ggplot(prdc.spp, aes(x = MMM_dev, y = Estimate, col = Reef)) +
  geom_line() +
  facet_grid(catBleaching~Species) +
  theme_light()

ggplot(prdc.spp, aes(x = MMM_dev, y = Estimate, col = catBleaching)) +
  geom_line() +
  facet_grid(Reef~Species) +
  theme_light()

d.Species:catBleaching - The interaction is significant e.Species:Reef - The interaction is significant

Species-specific models

Due to some species not being present for some reefs, we will run the species separately. The difficulty now is that not all bleaching categories will be present for each reef in each species. We can compare the results obtained from the 3-way interaction with the results obtained for each individuals species modelled separately.

First, we’ll have to subset the data to remove reefs with small sample sizes for each reef

#How many samples from each Species/Reef
psii %>% group_by(Species, Reef) %>% 
  dplyr::summarise(n = dplyr::n()) %>% 
  spread(key = Species, value = n)
## `summarise()` has grouped output by 'Species'. You can override using the
## `.groups` argument.
## # A tibble: 9 × 4
##   Reef         `A. cf humilis` `P. meandrina` `P. verrucosa`
##   <fct>                  <int>          <int>          <int>
## 1 Bougainville              73             36             59
## 2 Chilcott                  70             39             11
## 3 Flinders                 112             20             97
## 4 Frederick                 62             48              6
## 5 Herald                    62             16             36
## 6 Lihou                    117             23             71
## 7 Moore                     67             63             37
## 8 Saumarez                  NA             85             NA
## 9 Wreck                     96             52             24
ahum.1 = psii %>% filter(Species == "A. cf humilis") %>% 
  as.data.frame() %>% 
  droplevels()

pmea.1 = psii %>% filter(Species == "P. meandrina") %>% 
  droplevels()

pver.1 = psii %>% 
  filter(Species == "P. verrucosa") %>% 
  filter(! Reef %in% c("Frederick", "Chilcott")) %>% 
  droplevels()

What parameters to estimate?

Not all bleaching categories present at each reef, explore model selection with a subset of reefs

ahum.1 %>% group_by(Reef, catBleaching) %>% 
  dplyr::summarise(n = dplyr::n()) %>% 
  spread(key = catBleaching, value = n)
## `summarise()` has grouped output by 'Reef'. You can override using the
## `.groups` argument.
## # A tibble: 8 × 4
## # Groups:   Reef [8]
##   Reef         Healthy  Fair  Poor
##   <fct>          <int> <int> <int>
## 1 Bougainville       4    34    35
## 2 Chilcott          24    40     6
## 3 Flinders          38    66     8
## 4 Frederick         25    37    NA
## 5 Herald            18    35     9
## 6 Lihou             27    72    18
## 7 Moore             11    46    10
## 8 Wreck             52    44    NA
ahum.2 = ahum.1 %>% filter(! Reef %in% c("Frederick", "Wreck")) %>% 
  as.data.frame() %>% droplevels()
#model d and e? 
# model with 
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = ahum.2,
            curveid = Reef,
            pmodels = c( ~ 1,  ~ 1,  ~ Reef), bcVal = 0.5)

reef.ahum.1 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ 1, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = ahum.2,
            curveid = Reef,
            pmodels = c( ~ 1,  ~ Reef,  ~ Reef), bcVal = 0.5)

reef.ahum.2 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = ahum.2,
            curveid = Reef,
            pmodels = c( ~ Reef,  ~ Reef,  ~ Reef), bcVal = 0.5)
#reef.ahum.3 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
#                    random = e ~ 1|Vial,
#                    fixed = list(b ~ Reef, d ~ Reef, e ~ Reef),
#                    start = coef(modNaive.ahum), 
#                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(reef.ahum.1, reef.ahum.2)
##             df       AIC
## reef.ahum.1 10 -1177.798
## reef.ahum.2 15 -1221.154

Can’t run any models to estimate b, we can assume that the gradient is the same within species. However, we get marginal improvements with d and e. Best model: reef.ahum.2

Accounting for bleaching

How does the catBleaching affect e and d

modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = ahum.2,
            curveid = Reef:catBleaching,
            pmodels = c( ~ 1,  ~ Reef*catBleaching,  ~ Reef*catBleaching), bcVal = 0.5)

reef.ahum.3 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ Reef*catBleaching, e ~ Reef*catBleaching),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = ahum.2,
            curveid = Reef:catBleaching,
            pmodels = c( ~ 1,  ~ Reef,  ~ Reef*catBleaching), bcVal = 0.5)

reef.ahum.4 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ Reef, e ~ Reef*catBleaching),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = ahum.2,
            curveid = Reef:catBleaching,
            pmodels = c( ~ 1,  ~ Reef*catBleaching,  ~ Reef), bcVal = 0.5)

reef.ahum.5 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ Reef*catBleaching, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))


AIC(reef.ahum.2, reef.ahum.3, reef.ahum.4, reef.ahum.5)
##             df       AIC
## reef.ahum.2 15 -1221.154
## reef.ahum.3 39 -1220.852
## reef.ahum.4 27 -1221.605
## reef.ahum.5 27 -1225.515
anova(reef.ahum.2, reef.ahum.5)
##             Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## reef.ahum.2     1 15 -1221.154 -1157.905 625.5770                        
## reef.ahum.5     2 27 -1225.515 -1111.667 639.7575 1 vs 2 28.36111  0.0049
summary(reef.ahum.5)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: median.yield ~ NLS.L3(MMM_dev, b, d, e) 
##  Data: ahum.2 
##         AIC       BIC   logLik
##   -1225.515 -1111.667 639.7575
## 
## Random effects:
##  Formula: e ~ 1 | Vial
##         e.(Intercept)   Residual
## StdDev:  1.864301e-05 0.06748164
## 
## Fixed effects: list(b ~ 1, d ~ Reef * catBleaching, e ~ Reef) 
##                                     Value  Std.Error  DF   t-value p-value
## b                               -1.351935 0.05987450 336 -22.57947  0.0000
## d.(Intercept)                    0.544895 0.04303904 336  12.66048  0.0000
## d.ReefChilcott                   0.091942 0.04669440 336   1.96901  0.0498
## d.ReefFlinders                   0.055393 0.04537504 336   1.22077  0.2230
## d.ReefHerald                     0.077505 0.04804852 336   1.61305  0.1077
## d.ReefLihou                      0.081456 0.04623146 336   1.76192  0.0790
## d.ReefMoore                      0.068564 0.05027447 336   1.36378  0.1735
## d.catBleachingFair               0.027588 0.04510576 336   0.61162  0.5412
## d.catBleachingPoor              -0.042135 0.04500108 336  -0.93630  0.3498
## d.ReefChilcott:catBleachingFair -0.035716 0.05022179 336  -0.71116  0.4775
## d.ReefFlinders:catBleachingFair -0.018777 0.04833243 336  -0.38850  0.6979
## d.ReefHerald:catBleachingFair   -0.056621 0.05192923 336  -1.09036  0.2763
## d.ReefLihou:catBleachingFair    -0.043119 0.04903014 336  -0.87944  0.3798
## d.ReefMoore:catBleachingFair    -0.068826 0.05326301 336  -1.29219  0.1972
## d.ReefChilcott:catBleachingPoor  0.022370 0.05755472 336   0.38868  0.6978
## d.ReefFlinders:catBleachingPoor  0.043943 0.05610708 336   0.78320  0.4341
## d.ReefHerald:catBleachingPoor    0.002716 0.05683793 336   0.04778  0.9619
## d.ReefLihou:catBleachingPoor     0.028739 0.05197094 336   0.55298  0.5806
## d.ReefMoore:catBleachingPoor    -0.077685 0.05687225 336  -1.36595  0.1729
## e.(Intercept)                    6.917471 0.10509048 336  65.82395  0.0000
## e.ReefChilcott                  -0.091070 0.14820155 336  -0.61450  0.5393
## e.ReefFlinders                   0.040367 0.14753203 336   0.27362  0.7845
## e.ReefHerald                    -0.536771 0.15546134 336  -3.45276  0.0006
## e.ReefLihou                      0.057506 0.14870335 336   0.38671  0.6992
## e.ReefMoore                     -0.144701 0.16376371 336  -0.88359  0.3775
##  Correlation: 
##                                 b      d.(In) d.RfCh d.RfFl d.RfHr d.RfLh
## d.(Intercept)                    0.046                                   
## d.ReefChilcott                  -0.012 -0.920                            
## d.ReefFlinders                  -0.022 -0.948  0.873                     
## d.ReefHerald                    -0.018 -0.895  0.824  0.848              
## d.ReefLihou                     -0.018 -0.930  0.856  0.882  0.832       
## d.ReefMoore                     -0.014 -0.855  0.788  0.811  0.765  0.796
## d.catBleachingFair              -0.006 -0.940  0.866  0.891  0.841  0.874
## d.catBleachingPoor              -0.017 -0.946  0.871  0.897  0.847  0.880
## d.ReefChilcott:catBleachingFair  0.000  0.844 -0.910 -0.800 -0.756 -0.785
## d.ReefFlinders:catBleachingFair  0.003  0.877 -0.808 -0.920 -0.785 -0.816
## d.ReefHerald:catBleachingFair    0.003  0.816 -0.752 -0.774 -0.908 -0.760
## d.ReefLihou:catBleachingFair     0.006  0.864 -0.796 -0.820 -0.774 -0.924
## d.ReefMoore:catBleachingFair     0.000  0.795 -0.733 -0.754 -0.712 -0.740
## d.ReefChilcott:catBleachingPoor -0.006  0.738 -0.803 -0.701 -0.662 -0.688
## d.ReefFlinders:catBleachingPoor  0.015  0.758 -0.699 -0.793 -0.679 -0.706
## d.ReefHerald:catBleachingPoor    0.010  0.749 -0.690 -0.710 -0.834 -0.697
## d.ReefLihou:catBleachingPoor     0.010  0.819 -0.754 -0.776 -0.733 -0.876
## d.ReefMoore:catBleachingPoor    -0.003  0.747 -0.689 -0.709 -0.670 -0.696
## e.(Intercept)                   -0.108 -0.126  0.112  0.117  0.110  0.114
## e.ReefChilcott                   0.064  0.088 -0.154 -0.083 -0.078 -0.081
## e.ReefFlinders                   0.205  0.095 -0.082 -0.150 -0.081 -0.084
## e.ReefHerald                     0.047  0.084 -0.076 -0.078 -0.144 -0.077
## e.ReefLihou                      0.192  0.094 -0.081 -0.085 -0.080 -0.145
## e.ReefMoore                      0.131  0.083 -0.073 -0.076 -0.072 -0.074
##                                 d.RfMr d.ctBF d.ctBP d.RC:BF d.RF:BF d.RH:BF
## d.(Intercept)                                                               
## d.ReefChilcott                                                              
## d.ReefFlinders                                                              
## d.ReefHerald                                                                
## d.ReefLihou                                                                 
## d.ReefMoore                                                                 
## d.catBleachingFair               0.804                                      
## d.catBleachingPoor               0.809  0.898                               
## d.ReefChilcott:catBleachingFair -0.722 -0.898 -0.806                        
## d.ReefFlinders:catBleachingFair -0.750 -0.933 -0.838  0.838                 
## d.ReefHerald:catBleachingFair   -0.698 -0.869 -0.780  0.780   0.811         
## d.ReefLihou:catBleachingFair    -0.740 -0.920 -0.826  0.826   0.859   0.799 
## d.ReefMoore:catBleachingFair    -0.921 -0.847 -0.760  0.761   0.790   0.736 
## d.ReefChilcott:catBleachingPoor -0.632 -0.702 -0.782  0.738   0.655   0.610 
## d.ReefFlinders:catBleachingPoor -0.649 -0.720 -0.802  0.647   0.742   0.625 
## d.ReefHerald:catBleachingPoor   -0.641 -0.711 -0.792  0.638   0.663   0.766 
## d.ReefLihou:catBleachingPoor    -0.700 -0.777 -0.866  0.698   0.725   0.675 
## d.ReefMoore:catBleachingPoor    -0.876 -0.710 -0.791  0.638   0.663   0.617 
## e.(Intercept)                    0.105  0.010  0.039 -0.008  -0.009  -0.008 
## e.ReefChilcott                  -0.074 -0.007 -0.027  0.017   0.006   0.006 
## e.ReefFlinders                  -0.076 -0.008 -0.030  0.006   0.015   0.006 
## e.ReefHerald                    -0.070 -0.007 -0.026  0.006   0.006   0.016 
## e.ReefLihou                     -0.076 -0.008 -0.029  0.006   0.007   0.006 
## e.ReefMoore                     -0.166 -0.007 -0.026  0.005   0.006   0.006 
##                                 d.RL:BF d.RM:BF d.RC:BP d.RF:BP d.RH:BP d.RL:BP
## d.(Intercept)                                                                  
## d.ReefChilcott                                                                 
## d.ReefFlinders                                                                 
## d.ReefHerald                                                                   
## d.ReefLihou                                                                    
## d.ReefMoore                                                                    
## d.catBleachingFair                                                             
## d.catBleachingPoor                                                             
## d.ReefChilcott:catBleachingFair                                                
## d.ReefFlinders:catBleachingFair                                                
## d.ReefHerald:catBleachingFair                                                  
## d.ReefLihou:catBleachingFair                                                   
## d.ReefMoore:catBleachingFair     0.779                                         
## d.ReefChilcott:catBleachingPoor  0.646   0.594                                 
## d.ReefFlinders:catBleachingPoor  0.662   0.610   0.627                         
## d.ReefHerald:catBleachingPoor    0.654   0.602   0.619   0.635                 
## d.ReefLihou:catBleachingPoor     0.821   0.658   0.677   0.694   0.685         
## d.ReefMoore:catBleachingPoor     0.653   0.814   0.619   0.634   0.626   0.685 
## e.(Intercept)                   -0.009  -0.008  -0.028  -0.031  -0.030  -0.033 
## e.ReefChilcott                   0.006   0.006   0.074   0.022   0.021   0.023 
## e.ReefFlinders                   0.007   0.006   0.019   0.020   0.023   0.025 
## e.ReefHerald                     0.006   0.005   0.019   0.021   0.044   0.022 
## e.ReefLihou                      0.007   0.006   0.019   0.024   0.023   0.038 
## e.ReefMoore                      0.006   0.021   0.018   0.021   0.020   0.022 
##                                 d.RM:BP e.(In) e.RfCh e.RfFl e.RfHr e.RfLh
## d.(Intercept)                                                             
## d.ReefChilcott                                                            
## d.ReefFlinders                                                            
## d.ReefHerald                                                              
## d.ReefLihou                                                               
## d.ReefMoore                                                               
## d.catBleachingFair                                                        
## d.catBleachingPoor                                                        
## d.ReefChilcott:catBleachingFair                                           
## d.ReefFlinders:catBleachingFair                                           
## d.ReefHerald:catBleachingFair                                             
## d.ReefLihou:catBleachingFair                                              
## d.ReefMoore:catBleachingFair                                              
## d.ReefChilcott:catBleachingPoor                                           
## d.ReefFlinders:catBleachingPoor                                           
## d.ReefHerald:catBleachingPoor                                             
## d.ReefLihou:catBleachingPoor                                              
## d.ReefMoore:catBleachingPoor                                              
## e.(Intercept)                   -0.029                                    
## e.ReefChilcott                   0.020  -0.708                            
## e.ReefFlinders                   0.020  -0.726  0.512                     
## e.ReefHerald                     0.020  -0.673  0.477  0.486              
## e.ReefLihou                      0.020  -0.719  0.508  0.537  0.481       
## e.ReefMoore                      0.100  -0.648  0.458  0.479  0.435  0.473
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.70090770 -0.37138277  0.01415396  0.46630647  7.08204470 
## 
## Number of Observations: 501
## Number of Groups: 141

There is a small effect of catBleaching on d but the interaction is not significant. We can see how catBleaching is affecting a subset of reefs, notably Bougainville and Moore. Bothe have small sample sizes of ‘Poor’ categories, which may explain the weak interaction. Overall, we can ignore the interaction and proceed with a simpler model that does not involve an interaction catBleaching*Reef. However, we may wish to control fo catBleaching as a random effect. Best model: reef.ahum.2

emmeans(reef.ahum.5, ~catBleaching|Reef, param = "d") %>% plot()

Furthermore, if we were to include the interaction, we have to remove some reefs where some bleaching categories are not present. We will favour including all reefs and incorporate catBleaching as a random effect.

modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = ahum.2,
            curveid = Reef,
            pmodels = c( ~ 1,  ~ Reef,  ~ Reef), bcVal = 0.5)

#random effects... can we improved the model

reef.ahum.6 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                   random = d ~ 1|catBleaching + e ~ 1|Vial,
                   fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
                   start = coef(modNaive.ahum), 
                   control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

reef.ahum.7 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = e ~ 1|catBleaching,
                    fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

reef.ahum.8 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = d ~ 1|catBleaching,
                    fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

reef.ahum.9 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = d + e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

reef.ahum.10 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = d + e ~ 1|catBleaching/Vial,
                    fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

reef.ahum.11 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
                    random = d + e ~ 1|catBleaching,
                    fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

AIC(reef.ahum.2, reef.ahum.6, reef.ahum.7, reef.ahum.8, reef.ahum.9, reef.ahum.10, reef.ahum.11)
##              df       AIC
## reef.ahum.2  15 -1221.154
## reef.ahum.6  15 -1221.154
## reef.ahum.7  15 -1221.154
## reef.ahum.8  15 -1221.154
## reef.ahum.9  17 -1217.154
## reef.ahum.10 20 -1211.154
## reef.ahum.11 17 -1217.154

Based on what we saw in the species model, interaction with catBleaching above and random effects here, catBleaching may be having most of an impact on d but not e but we still want to control for Vial on e so reef.ahum.6 would be the top choice. Best: reef.ahum.6

However, the models don’t converge when consider all reefs tested with each species. So we have to simply the model further. For consistency and to compare results, we will also need to use the same model for every species.

A. cf humilis

#with all reefs 
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = ahum.1,
            curveid = Reef,
            pmodels = c( ~ 1,  ~ 1,  ~ Reef), bcVal = 0.5)

reef.ahum = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.1,
                    random = d ~ 1|catBleaching + e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ 1, e ~ Reef),
                    start = coef(modNaive.ahum), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))

To protect our original models, we’ve saved them as a Rdata file

load("3_outputs/nlme_A.humilis by reef.RData")
emmeans(reef.ahum, ~Reef, param = "e") %>% plot(comparisons = TRUE)

emmeans(reef.ahum, pairwise~Reef, param = "e")$contrast %>% plot()

emmeans(reef.ahum, ~Reef, param = 'e', level = 0.95) 
##  Reef         emmean     SE  df lower.CL upper.CL
##  Bougainville   6.70 0.0941 468     6.52     6.89
##  Chilcott       6.91 0.1107 468     6.70     7.13
##  Flinders       6.93 0.0986 468     6.74     7.13
##  Frederick      6.80 0.1194 468     6.56     7.03
##  Herald         6.37 0.1156 468     6.14     6.59
##  Lihou          7.01 0.1050 468     6.80     7.22
##  Moore          6.64 0.1112 468     6.42     6.86
##  Wreck          8.26 0.0885 468     8.08     8.43
## 
## Confidence level used: 0.95
emmeans(reef.ahum, pairwise~Reef, param = "e")$contrast
##  contrast                 estimate    SE  df t.ratio p.value
##  Bougainville - Chilcott   -0.2082 0.144 468  -1.447 0.8348 
##  Bougainville - Flinders   -0.2268 0.135 468  -1.677 0.7024 
##  Bougainville - Frederick  -0.0941 0.151 468  -0.625 0.9985 
##  Bougainville - Herald      0.3386 0.148 468   2.288 0.3021 
##  Bougainville - Lihou      -0.3054 0.140 468  -2.186 0.3624 
##  Bougainville - Moore       0.0617 0.145 468   0.426 0.9999 
##  Bougainville - Wreck      -1.5534 0.127 468 -12.211 <.0001 
##  Chilcott - Flinders       -0.0186 0.147 468  -0.127 1.0000 
##  Chilcott - Frederick       0.1141 0.161 468   0.708 0.9968 
##  Chilcott - Herald          0.5468 0.159 468   3.441 0.0145 
##  Chilcott - Lihou          -0.0972 0.151 468  -0.644 0.9982 
##  Chilcott - Moore           0.2699 0.156 468   1.733 0.6654 
##  Chilcott - Wreck          -1.3452 0.140 468  -9.594 <.0001 
##  Flinders - Frederick       0.1327 0.153 468   0.869 0.9886 
##  Flinders - Herald          0.5655 0.151 468   3.746 0.0049 
##  Flinders - Lihou          -0.0786 0.139 468  -0.565 0.9992 
##  Flinders - Moore           0.2885 0.145 468   1.984 0.4942 
##  Flinders - Wreck          -1.3266 0.135 468  -9.853 <.0001 
##  Frederick - Herald         0.4328 0.165 468   2.624 0.1499 
##  Frederick - Lihou         -0.2113 0.157 468  -1.349 0.8794 
##  Frederick - Moore          0.1558 0.162 468   0.965 0.9791 
##  Frederick - Wreck         -1.4592 0.147 468  -9.899 <.0001 
##  Herald - Lihou            -0.6440 0.155 468  -4.155 0.0010 
##  Herald - Moore            -0.2769 0.160 468  -1.735 0.6641 
##  Herald - Wreck            -1.8920 0.144 468 -13.117 <.0001 
##  Lihou - Moore              0.3671 0.150 468   2.451 0.2194 
##  Lihou - Wreck             -1.2480 0.139 468  -8.982 <.0001 
##  Moore - Wreck             -1.6151 0.143 468 -11.299 <.0001 
## 
## P value adjustment: tukey method for comparing a family of 8 estimates
#Some significant differences

P. meandrina

#with all reefs 
modNaive.pmea <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = pmea.1,
            curveid = Reef,
            pmodels = c( ~ 1,  ~ 1,  ~ Reef), bcVal = 0.5)


reef.pmea = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = pmea.1,
                    random = d ~ 1|catBleaching + e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ 1, e ~ Reef),
                    start = coef(modNaive.pmea), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
load("3_outputs/nlme_P.meandrina by reef.RData")
emmeans(reef.pmea, ~Reef, param = "e") %>% plot(comparisons = TRUE)

emmeans(reef.pmea, pairwise~Reef, param = "e")$contrast %>% plot()

emmeans(reef.pmea, ~Reef, param = 'e', level = 0.95) 
##  Reef         emmean     SE  df lower.CL upper.CL
##  Bougainville   7.09 0.0973 271     6.89     7.28
##  Chilcott       7.08 0.1148 271     6.85     7.31
##  Flinders       8.11 0.1555 271     7.81     8.42
##  Frederick      7.06 0.1166 271     6.83     7.29
##  Herald         7.25 0.1885 271     6.88     7.62
##  Lihou          6.96 0.1770 271     6.61     7.31
##  Moore          7.26 0.1159 271     7.03     7.49
##  Saumarez       7.86 0.1167 271     7.63     8.09
##  Wreck          8.03 0.1016 271     7.83     8.23
## 
## Confidence level used: 0.95
emmeans(reef.pmea, pairwise~Reef, param = "e")$contrast
##  contrast                 estimate    SE  df t.ratio p.value
##  Bougainville - Chilcott   0.00496 0.150 271  0.033  1.0000 
##  Bougainville - Flinders  -1.02701 0.180 271 -5.717  <.0001 
##  Bougainville - Frederick  0.02894 0.148 271  0.195  1.0000 
##  Bougainville - Herald    -0.16248 0.211 271 -0.771  0.9975 
##  Bougainville - Lihou      0.12858 0.203 271  0.633  0.9994 
##  Bougainville - Moore     -0.17636 0.149 271 -1.184  0.9593 
##  Bougainville - Saumarez  -0.77570 0.146 271 -5.298  <.0001 
##  Bougainville - Wreck     -0.94347 0.136 271 -6.951  <.0001 
##  Chilcott - Flinders      -1.03197 0.194 271 -5.318  <.0001 
##  Chilcott - Frederick      0.02398 0.163 271  0.147  1.0000 
##  Chilcott - Herald        -0.16743 0.220 271 -0.762  0.9977 
##  Chilcott - Lihou          0.12362 0.208 271  0.594  0.9996 
##  Chilcott - Moore         -0.18131 0.161 271 -1.123  0.9704 
##  Chilcott - Saumarez      -0.78065 0.164 271 -4.758  0.0001 
##  Chilcott - Wreck         -0.94842 0.154 271 -6.144  <.0001 
##  Flinders - Frederick      1.05595 0.188 271  5.603  <.0001 
##  Flinders - Herald         0.86454 0.243 271  3.562  0.0128 
##  Flinders - Lihou          1.15559 0.241 271  4.792  0.0001 
##  Flinders - Moore          0.85065 0.191 271  4.445  0.0004 
##  Flinders - Saumarez       0.25131 0.184 271  1.366  0.9095 
##  Flinders - Wreck          0.08355 0.176 271  0.475  0.9999 
##  Frederick - Herald       -0.19141 0.220 271 -0.871  0.9943 
##  Frederick - Lihou         0.09964 0.214 271  0.465  0.9999 
##  Frederick - Moore        -0.20529 0.161 271 -1.273  0.9382 
##  Frederick - Saumarez     -0.80463 0.157 271 -5.132  <.0001 
##  Frederick - Wreck        -0.97240 0.147 271 -6.605  <.0001 
##  Herald - Lihou            0.29105 0.258 271  1.128  0.9695 
##  Herald - Moore           -0.01388 0.219 271 -0.063  1.0000 
##  Herald - Saumarez        -0.61322 0.219 271 -2.797  0.1212 
##  Herald - Wreck           -0.78099 0.212 271 -3.678  0.0085 
##  Lihou - Moore            -0.30493 0.211 271 -1.446  0.8788 
##  Lihou - Saumarez         -0.90427 0.218 271 -4.147  0.0015 
##  Lihou - Wreck            -1.07204 0.211 271 -5.090  <.0001 
##  Moore - Saumarez         -0.59934 0.160 271 -3.737  0.0069 
##  Moore - Wreck            -0.76711 0.151 271 -5.082  <.0001 
##  Saumarez - Wreck         -0.16777 0.141 271 -1.187  0.9587 
## 
## P value adjustment: tukey method for comparing a family of 9 estimates
#Some significant differences

P. verrucosa

#with all reefs 
modNaive.pver <- drm(median.yield ~ MMM_dev, fct = L.3(), 
            data = pver.1,
            curveid = Reef,
            pmodels = c( ~ 1,  ~ 1,  ~ Reef), bcVal = 0.5)


reef.pver = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = pver.1,
                    random = d ~ 1|catBleaching + e ~ 1|Vial,
                    fixed = list(b ~ 1, d ~ 1, e ~ Reef),
                    start = coef(modNaive.pver), 
                    control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
load("3_outputs/nlme_P.verrucosa by reef.RData")
emmeans(reef.pver, ~Reef, param = "e") %>% plot(comparisons = TRUE)

emmeans(reef.pver, pairwise~Reef, param = "e")$contrast %>% plot()

emmeans(reef.pver, ~Reef, param = 'e', level = 0.95) 
##  Reef         emmean     SE  df lower.CL upper.CL
##  Bougainville   7.51 0.0587 229     7.40     7.63
##  Flinders       8.11 0.0903 229     7.93     8.28
##  Herald         7.26 0.1251 229     7.01     7.51
##  Lihou          7.74 0.1307 229     7.48     7.99
##  Moore          7.54 0.1373 229     7.27     7.81
##  Wreck          8.09 0.1279 229     7.83     8.34
## 
## Confidence level used: 0.95
emmeans(reef.pver, pairwise~Reef, param = "e")$contrast
##  contrast                estimate     SE  df t.ratio p.value
##  Bougainville - Flinders  -0.5949 0.0938 229 -6.339  <.0001 
##  Bougainville - Herald     0.2514 0.1373 229  1.830  0.4483 
##  Bougainville - Lihou     -0.2240 0.1311 229 -1.709  0.5275 
##  Bougainville - Moore     -0.0317 0.1418 229 -0.224  0.9999 
##  Bougainville - Wreck     -0.5736 0.1311 229 -4.375  0.0003 
##  Flinders - Herald         0.8463 0.1536 229  5.511  <.0001 
##  Flinders - Lihou          0.3708 0.1314 229  2.822  0.0576 
##  Flinders - Moore          0.5631 0.1481 229  3.801  0.0025 
##  Flinders - Wreck          0.0212 0.1348 229  0.157  1.0000 
##  Herald - Lihou           -0.4754 0.1798 229 -2.644  0.0910 
##  Herald - Moore           -0.2832 0.1844 229 -1.535  0.6418 
##  Herald - Wreck           -0.8250 0.1785 229 -4.623  0.0001 
##  Lihou - Moore             0.1923 0.1728 229  1.112  0.8759 
##  Lihou - Wreck            -0.3496 0.1611 229 -2.170  0.2559 
##  Moore - Wreck            -0.5419 0.1746 229 -3.104  0.0259 
## 
## P value adjustment: tukey method for comparing a family of 6 estimates
#Some significant differences

Absolute Temperature ED50s (Table S5 SOM)

ahum.abs = ahum.1 %>% group_by(Reef) %>% 
  summarise(MMM = mean(MMM)) %>% 
  left_join(emmeans(reef.ahum, ~Reef, param = "e") %>% as.data.frame(), by = "Reef") %>%  
  mutate(absolute = MMM + emmean) %>% 
  mutate(absolute.lowerCL = MMM + lower.CL) %>% 
  mutate(absolute.upperCL = MMM + upper.CL)

pmea.abs = pmea.1%>% group_by(Reef) %>% 
  summarise(MMM = mean(MMM)) %>% 
  left_join(emmeans(reef.pmea, ~Reef, param = "e") %>% as.data.frame(), by = "Reef") %>%  
  mutate(absolute = MMM + emmean) %>% 
  mutate(absolute.lowerCL = MMM + lower.CL) %>% 
  mutate(absolute.upperCL = MMM + upper.CL)

pver.abs = pver.1 %>% group_by(Reef) %>% 
  summarise(MMM = mean(MMM)) %>% 
  left_join(emmeans(reef.pver, ~Reef, param = "e") %>% as.data.frame(), by = "Reef") %>%  
  mutate(absolute = MMM + emmean) %>% 
  mutate(absolute.lowerCL = MMM + lower.CL) %>% 
  mutate(absolute.upperCL = MMM + upper.CL)

2.5 Model predictions: Species*Reef

##Create predictions data frame
predframe.ahum = with(ahum.1, expand.grid(Reef = levels(Reef), MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))

#Adding confidence intervals
prdc.ahum = as.data.frame(nlraa::predict_nlme(reef.ahum, newdata = predframe.ahum, interval = "confidence", level = 0.95))
prdc.ahum = bind_cols(predframe.ahum, prdc.ahum)

# Predict ED50 for each species 
#(Since we're not varyin d, we don't have to do this but could be handy nevertheless. Estimate values should not change much between reefs)
ED50.ahum = emmeans(reef.ahum, ~Reef, param= "e") %>% as_tibble() %>% 
  dplyr::select(Reef, emmean) %>% 
  dplyr::rename(MMM_dev = emmean)
  
prdc.ED50.ahum = as.data.frame(nlraa::predict_nlme(reef.ahum, newdata = ED50.ahum, interval = "confidence", level = 0.95))
prdc.ED50.ahum = bind_cols(ED50.ahum, prdc.ED50.ahum)

P. meandrina

##Create predictions data frame
predframe.pmea = with(pmea.1, expand.grid(Reef = levels(Reef), MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))

#Adding confidence intervals
prdc.pmea = as.data.frame(nlraa::predict_nlme(reef.pmea, newdata = predframe.pmea, interval = "confidence", level = 0.95))
prdc.pmea = bind_cols(predframe.pmea, prdc.pmea)


# Predict ED50 for each species
ED50.pmea = emmeans(reef.pmea, ~Reef, param= "e") %>% as_tibble() %>% 
  dplyr::select(Reef, emmean) %>% 
  dplyr::rename(MMM_dev = emmean)
  
prdc.ED50.pmea = as.data.frame(nlraa::predict_nlme(reef.pmea, newdata = ED50.pmea, interval = "confidence", level = 0.95))
prdc.ED50.pmea = bind_cols(ED50.pmea, prdc.ED50.pmea)

P. verrucosa

##Create predictions data frame
predframe.pver = with(pver.1, 
                     expand.grid(Reef = levels(Reef),
                                 MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))

#Adding confidence intervals
prdc.pver = as.data.frame(nlraa::predict_nlme(reef.pver, newdata = predframe.pver, 
                                           interval = "confidence", level = 0.95))
prdc.pver = bind_cols(predframe.pver, prdc.pver)


# Predict ED50 for each species
ED50.pver = emmeans(reef.pver, ~Reef, param= "e") %>% as_tibble() %>% 
  dplyr::select(Reef, emmean) %>% 
  dplyr::rename(MMM_dev = emmean)
  
prdc.ED50.pver = as.data.frame(nlraa::predict_nlme(reef.pver, newdata = ED50.pver, interval = "confidence", level = 0.95))
prdc.ED50.pver = bind_cols(ED50.pver, prdc.ED50.pver)

2.6 Figure: Reef comparison

#Base order of reefs and colours on the average ED50.
Reef.ED50 = prdc.ED50.ahum %>% mutate(Species = "A. cf humilis") %>% dplyr::select(Species, Reef, MMM_dev) %>% 
  bind_rows(prdc.ED50.pmea %>% mutate(Species = "P. meandrina") %>% dplyr::select(Species, Reef, MMM_dev)) %>% 
  bind_rows(prdc.ED50.pver %>% mutate(Species = "P. verrucosa") %>% dplyr::select(Species, Reef, MMM_dev)) %>% 
  spread(key = Species, value = MMM_dev) %>% 
  rowwise(Reef) %>% 
  dplyr::summarise(average = mean(c(`A. cf humilis`, `P. meandrina`, `P. verrucosa`), na.rm = T)) %>% 
  arrange(average) 
## `summarise()` has grouped output by 'Reef'. You can override using the
## `.groups` argument.
Reef.order = fct_reorder(Reef.ED50$Reef, Reef.ED50$average)

#Combine all data and reorder Reef factor
line.pred.dat = prdc.ahum %>% mutate(Species = "A. cf humilis") %>% 
  bind_rows(prdc.pmea %>% mutate(Species = "P. meandrina")) %>% 
  bind_rows(prdc.pver %>% mutate(Species = "P. verrucosa")) %>% 
  mutate(Reef = factor(Reef, levels(Reef.order)))
  
point.obs.dat = ahum.1 %>% 
  bind_rows(pmea.1) %>% 
  bind_rows(pver.1) %>%
  mutate(Reef = factor(Reef, levels(Reef.order)))

point.ED50.dat = prdc.ED50.ahum %>% mutate(Species = "A. cf humilis") %>% 
  bind_rows(prdc.ED50.pmea %>% mutate(Species = "P. meandrina")) %>% 
  bind_rows(prdc.ED50.pver %>% mutate(Species = "P. verrucosa")) %>% 
  mutate(Reef = factor(Reef, levels(Reef.order)))

Fig 2B: PSII comparison by Reef

col = wes_palette("Zissou1", 9, type = "continuous")

Reef.plot = ggplot() +
  geom_segment(data = point.ED50.dat, aes(x = MMM_dev, xend = MMM_dev, y = -Inf, yend = Estimate , color = Reef), 
                  linetype = 1, size = .3) +
  geom_line(data = line.pred.dat, aes(colour=Reef, y=Estimate, x=MMM_dev), size = .3) +
  geom_point(data = point.ED50.dat, aes(x = MMM_dev, y = Estimate, fill = Reef), size = 1.5, stroke = .2, shape = 21, col = "white") +
  geom_text(data = data.frame(Species = c("A. cf humilis", "P. meandrina", "P. verrucosa")),
            aes(y = .71, x = 2, label = Species), size = 3, fontface = "italic", family = "Helvetica", hjust = 0) +
  #scale_y_continuous(labels = scales::percent, breaks = c(0, .5, 1)) +
  scale_x_continuous(breaks = c(3,6,9)) +
  coord_cartesian(xlim = c(2,9.6), ylim = c(0,.75)) +
  scale_fill_manual(values=col) + 
  scale_colour_manual(values=col) +
  facet_grid(Species~.) +
  labs(x="Temp. above local MMM (°C) ", 
       y="PSII-yield") +
  theme_classic() +
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        strip.text = element_blank(),
        strip.background = element_blank(),
        legend.title = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        legend.position = "right", 
        legend.text = element_text(size = 7, family = "Helvetica"),
        legend.key.size = unit(2, units = "mm"),
        legend.margin = margin(0,0,0,0, unit = "mm"),
        axis.text.y = element_text(size = 8, family = "Helvetica"),
        axis.title = element_text(size = 9, family = "Helvetica"))
  
Reef.plot

2.6 Model comparison

Does it make a difference to estimate e using the 3 way interaction of for each reef separately?

Fig SX: model comparison

#All species in the same model
mod.all = emmeans(reef.1, ~Reef|Species, param = "e")  %>% as.tibble() %>% 
  dplyr::select(Species, Reef, emmean, lower.CL, upper.CL)
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
#Species modelled separately
mod.sep = emmeans(reef.ahum, ~Reef, param= "e") %>% as.tibble() %>% 
              mutate(Species = "A. cf humilis") %>% dplyr::select(Species, Reef, emmean, lower.CL, upper.CL) %>% 
  bind_rows(emmeans(reef.pmea, ~Reef, param= "e") %>% as.tibble() %>% 
              mutate(Species = "P. meandrina") %>% dplyr::select(Species, Reef, emmean, lower.CL, upper.CL)) %>% 
  bind_rows(emmeans(reef.pver, ~Reef, param= "e") %>% as.tibble() %>% 
              mutate(Species = "P. verrucosa") %>% dplyr::select(Species, Reef, emmean, lower.CL, upper.CL)) %>% 
  mutate(Reef = factor(Reef, levels(Reef.order)))

ggplot(bind_rows(mod.all %>% mutate(model = "combined"), mod.sep %>% mutate(model = "separate")), 
       aes(y = emmean, x = Reef, ymin = lower.CL, ymax = upper.CL, col = model)) +
  geom_pointrange(position = position_dodge(width = .5)) +
  scale_color_manual(values = c("grey30", "skyblue2")) +
  facet_wrap(~Species, scales = "free_x") +
  labs(y = "ED50 estimate", x = "") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.line = element_blank(), 
        panel.border = element_rect(color = "black", fill = "transparent", size = 1),
        strip.background = element_blank(), 
        strip.text = element_text(face = "italic"),
        legend.title = element_blank(),
        legend.position = c(.9,.15),
        legend.background = element_rect(fill = "transparent"))

ggsave("2_figures/FigS4_Model comparison.png", width = 140, height = 70, units = "mm")

It makes no difference what model we use.

3. DRC Results

Table S1. ED50 values for each reef nested within species + 95 CI

library(gt)
library(tidyverse)

tableED50reef = mod.sep %>% 
  mutate(emmean = round(emmean,2), CI = paste(round(lower.CL,2), round(upper.CL,2), sep = "-")) %>% 
  dplyr::select(-lower.CL, -upper.CL) %>% 
  mutate(ED50 = paste(emmean, CI, sep = " ")) %>% 
  dplyr::select(-emmean, -CI) %>% 
  spread(key = Species, value = ED50) %>% 
  mutate(Reef = factor(Reef, levels = Reef.order)) %>% arrange(Reef)

gt_tbl.reefED50 <- gt(tableED50reef) %>% 
  tab_header(
  title = "Table S4. Summary tables derived from emmeans for ED50 values related to each reef and species combination (95% confidence intervals). ",) %>% 
  tab_spanner(label = "Species",
    columns = c("A. cf humilis", "P. meandrina", "P. verrucosa"))

#gt_tbl.reefED50

Figure 2

g1 = ggarrange(
  #Species plot
  Species.plot + labs(x = "") + 
    theme(axis.text.x = element_blank(), 
          plot.margin  = margin(0.1,0.1,0,.6, unit = "cm")), 
  #Reef plot        
  Reef.plot + theme(plot.margin  = margin(0.1,0.24,0.1,.6, unit = "cm")), 
          nrow = 2, heights = c(.33,.66), labels = c("a", "b"), 
          font.label = list(size = 12, family = "Helvetica")) 

g1

ggsave("2_figures/Fig2_v1.png", dpi =300, width=90, height=110, units = "mm")

#4. Drivers of heat tolerance among reefs in the CSMP

There are significant differences in the estimated ED50 among reefs. Can different climatic regimes be driving these patterns?

4.2 Data

Fig 3a . PSII-50 by reef and species with 95%CI

col = wes_palette("GrandBudapest1", 4, type = "continuous")

Fig3a =  ggplot(mod.sep, aes(y = emmean, x = Reef, ymin = lower.CL, ymax = upper.CL, col = Species)) +
  geom_pointrange(position = position_dodge(width = .4), size = .2) +
  scale_y_continuous(expand = c(0,0), lim = c(5.9,8.5), breaks = c(6.5, 7, 7.5, 8)) +
  scale_colour_manual(values=col) +
  labs(x="", 
       y="ED50 (95% CI)") +
  theme_classic() + 
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        legend.title = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        legend.position = c(.22,.9),
        legend.text = element_text(face = "italic", size = 10, family = "Helvetica"),
        legend.key.size = unit(3, units = "mm"),
        legend.margin = margin(0,0,0,0, unit = "mm"),
        axis.text = element_text(size = 14, family = "Helvetica"),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title = element_text(size = 14, family = "Helvetica"))

Fig3a

We can observe a large variation in the point of inflection between reefs in the CSMP that span several degrees of latitude and an array of climatic regimes.

Climate predictors of PSII-50

We will use the PSII-50 values derived from separate reef models and test against a suite of environmental variables

Thermal history metadata (based on NOAA ….. ) DHW3: Number of Degree Heating Weeks (DHW) events greater than 3 between 1985 - 2020 DHW4: Number of DHW events greater than 4C-weeks between 1985 - 2020 DHW6: Number of DHW events greater than 6C-weeks between 1985 - 2020 DHW8: Number of DHW events greater than 8C-weeks between 1985 - 2020 DHW9: Number of DHW events greater than 9C-weeks between 1985 - 2020 maxDHW: Maximum DHW experienced between 1985 - 2020 meanDHW: Average DHWs experienced between 1985 - 2020

MMM: Maximum Monthly Mean SST between 1986 - 2010 meanSST: Average Sea Surface Temperature (SST) between 1985 - 2020 maxSST: Maximum Sea Surface Temperature (SST) between 1985 - 2020 minSST: Max Sea Surface Temperature (SST) between 1985 - 2020 rangeSST: Range between Max and Min SST values between 1985 - 2020 varSST: Variability in SST between 1985 - 2020

Recent.maxDHW Maximum DHW exposue between 2016 - 2020 Recent.meanDHW Average maximum DHW exposure between 2016 - 2020

ReturnDHW3 The return time between events where DHW exceeded 3C-weeks ReturnDHW4 The return time between events where DHW exceeded 4C-weeks

Merging all data frames together

#Load data for thermal history data 
load("1_data/SiteDisturbanceHistory_DHW.RData")
load("1_data/reef.info.RData")

#DHW values for all reefs at the time of sampling (Feb-Mar 2020)
DHW2020 = bind_rows(ahum.1, pmea.1, pver.1) %>% 
  group_by(Reef, Site, DHW) %>% 
  dplyr::summarise() %>% 
  group_by(Reef) %>% 
  dplyr::summarise(DHW2020 = mean(DHW, na.rm = TRUE)) 
## `summarise()` has grouped output by 'Reef', 'Site'. You can override using the
## `.groups` argument.
#we merge data frames for thermal history metrics and geographic information with the ED50 values 
clim.variables = site.bleachings %>% ungroup() %>% dplyr::select(-Sector, -Site, -Month) %>% 
  group_by(Reef) %>% dplyr::summarise_all(mean) %>% 
  left_join(DHW2020, by = "Reef") %>% 
  left_join(reef.info %>% ungroup() %>% dplyr::select(-SectorRegion, -Sector, -Region), by = "Reef") %>% 
  filter(Reef %in% mod.sep$Reef)

psii.prediction = clim.variables %>% 
  gather(key = variable, value = value, -Reef) %>% 
  full_join(mod.sep %>% dplyr::rename(ED50 = emmean), by = "Reef") %>% 
  filter(! is.na(value)) %>% 
  dplyr::mutate(value = round(value, 2))

#linear relationships
col = wes_palette("GrandBudapest1", 4, type = "continuous")

ggplot(psii.prediction, aes(y = ED50, x = value, col = Species, fill = Species)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_fill_manual(values=col) + 
  scale_colour_manual(values=col) +
  labs(y = "ED50", x = "Environmental metric") +
  facet_wrap(~variable, scales = "free") +
  theme_classic() +
  theme(axis.line = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1),       
      axis.text = element_text(size = 8, family = "Helvetica"), 
      strip.text = element_text(size = 7)) 
## `geom_smooth()` using formula 'y ~ x'

ggsave("2_figures/FigS5_correlation plot.png", height = 5, width = 7, units = "in")
## `geom_smooth()` using formula 'y ~ x'

At first glance, we see that many variables show some relationship with ED50, but that heat tolerance drivers might be species-specific (ED50) in response to environmental variables. For example, DHW3 has strong correlation to ED50 value for A. humilis but the trend looks less strong for either of the Poci species, which may point towards an interaction.

4.1 Correlation matrix of thermal history metrics

Correlations between climate variables

Before we interpret these relationships, it’s important to consider that many of these variables are co-linear. Any metrics that have high correlation may represent the same climate variability and do not need to be included.

library(corrplot)
## corrplot 0.92 loaded
clim.variables %>% dplyr::select(-Reef) %>% cor() %>% corrplot()

Correlations between ED50 values and thermal variables

Let’s look at correlation coefficients from these data to see which variables show strong correlation.

psii.prediction %>% 
  group_by(variable, Species) %>% 
  summarise(correlation = cor(value, ED50)) %>% 
  ggplot(aes(y = variable, x = Species, fill = correlation)) +
  geom_tile() +
  scale_fill_gradient2(low = "darkred", mid = "white", high = "blue3") +
  theme_classic() + theme(axis.line = element_blank())

ggsave("2_figures/FigS7_correlation by species.png", height = 5, width = 7, units = "in")


save(psii.prediction, file = "1_data/psii.prediction.Rdata")

Looking at each species correlation values, we find that each species responds differently to each environmental variables. In common, all species seem to have a very fairly strong positive relationship to varSST, and negative relationship to recent meanDHW.

Here we are looking at the correlation values between enviro and ED50 selecting a subset with any variables that are > 0.40.

psii.prediction %>% 
  group_by(variable) %>% 
  summarise(correlation = cor(value, ED50))
## # A tibble: 24 × 2
##    variable   correlation
##    <chr>            <dbl>
##  1 Complexity     -0.202 
##  2 DHW2            0.353 
##  3 DHW2020        -0.501 
##  4 DHW3            0.353 
##  5 DHW4            0.458 
##  6 DHW6           -0.0234
##  7 DHW8           -0.375 
##  8 DHW9           -0.177 
##  9 Lat            -0.498 
## 10 Long            0.365 
## # … with 14 more rows
## # ℹ Use `print(n = ...)` to see more rows
psii.wide = psii.prediction %>% spread(variable, value)
save(psii.wide, file = "1_data/psii.wide.Rdata")

We can take a subset of variables to represent the range of climate regimes and that have low colinearity to other variables. We are excluding MMM and DHW2020 because the collinearity is > 0.80 to other variables that are of importance and have a stronger correlation to ED50s (e.g. recent.meanDHW).

climate.subset = c("DHW2020", "DHW4", "Lat", "meanSST", "minSST", "MMM", "rangeSST", "recent.meanDHW", "returnDHW6", "varSST")

clim.variables %>% dplyr::select(all_of(climate.subset)) %>% 
  cor()  %>% corrplot() 

#we will keep all but range SST. 

We are now subsetting variables for any that have high correlation to ED50. - First, we excluded latitude, meanSST, minSST, maxSST because they are all highly correlated to each other. - Then, we also exclude DHW2020 and MMM because they are highly correlated to recent.meanDHW. We choose to keep recent.meanDHW because it has the highest correlation with ED50 out of any of these 3 variables. - We then exclude rangeSST because of high correlation with varSST. VarSST has stronger relationship with ED50 so we preferentially keep it in. This leaves us with 4 remaining variables to go into the first dredge model: DHW4, recent.meanDHW, returnDHW6, and varSST.

thermal.subset = c("DHW4", "recent.meanDHW", "returnDHW6", "varSST")

clim.variables%>% dplyr::select(all_of(thermal.subset)) %>% 
  cor() %>% corrplot()

Below, the four plots represent some of the strongest climate variables that are not too strongly correlated with each other, but show some correlation with ED50. This allows us to assess whether some variables might have an interaction with species. We can use this subset to investigate the drivers of ED50 in a dredge model.

col = wes_palette("GrandBudapest1", 4, type = "continuous")

psii.prediction %>% filter(variable %in% thermal.subset) %>% 
  ggplot(aes(y = ED50, x = value, col = Species, fill = Species)) +
  geom_smooth(method = "lm") +
  geom_point() +
  scale_fill_manual(values=col) + 
  scale_colour_manual(values=col) +
  labs(y = "ED50") +
  facet_wrap(~variable, scales = "free") +
  theme_classic() +
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        axis.text = element_text(size = 8, family = "Helvetica"),
        axis.title = element_text(size = 9, family = "Helvetica"))
## `geom_smooth()` using formula 'y ~ x'

ggsave("2_figures/FigS8_EnviroRelationships.png", height = 5, width = 7, units = "in")
## `geom_smooth()` using formula 'y ~ x'

4.3 Dredging: Candidate Model Selection

Initial dredge

library(MuMIn)
library(lme4)
## Warning: package 'lme4' was built under R version 4.0.5
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
options(na.action = "na.fail")

#creating a data frame with only the 4 variables selected (DHW4, recent.meanDHW, returnDHW6, and varSST)
dredge.dat1 = psii.prediction %>% filter(variable %in% thermal.subset) %>% 
  spread(key = variable, value = value) %>% 
  dplyr::select(-ends_with(".CL"), -Reef)

#dredge model that accounts for species + the four variables in the model. 
dredge.mod1 <- dredge(lm(ED50~Species +., data = dredge.dat1), rank = "AIC", m.lim = c(1,4))
## Fixed term is "(Intercept)"
#just look at the top 10 scoring models 
head(dredge.mod1, 10)
## Global model call: lm(formula = ED50 ~ Species + ., data = dredge.dat1)
## ---
## Model selection table 
##    (Int)     DHW rcn.mDH     rDH Spc    vSS df logLik  AIC delta weight
## 16 6.607 0.12220 -0.1630 0.05856   +         7 -0.548 15.1  0.00  0.471
## 14 5.212 0.17540         0.08732   +         6 -2.592 17.2  2.09  0.166
## 30 4.335 0.14940         0.06693   + 0.4729  7 -2.277 18.6  3.46  0.084
## 12 7.998 0.06906 -0.2910           +         6 -3.664 19.3  4.23  0.057
## 28 5.601 0.06648 -0.1862           + 0.7482  7 -2.720 19.4  4.34  0.054
## 11 8.734         -0.3401           +         5 -4.907 19.8  4.72  0.044
## 27 6.183         -0.2278           + 0.7871  6 -3.967 19.9  4.84  0.042
## 15 8.330         -0.2959 0.02869   +         6 -4.125 20.2  5.15  0.036
## 26 2.637 0.08607                   + 1.4910  6 -4.335 20.7  5.57  0.029
## 31 6.697         -0.2356 0.01616   + 0.5583  7 -3.789 21.6  6.48  0.018
## Models ranked by AIC(x)
#extracting the best fit model using AIC score 
bestmodel <- get.models(dredge.mod1, 1)[[1]]
lm.dredge1 <- lm(bestmodel, data = dredge.dat1)
summary(lm.dredge1)
## 
## Call:
## lm(formula = bestmodel, data = dredge.dat1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44268 -0.15135  0.02352  0.12432  0.53704 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          6.60656    0.83588   7.904  4.3e-07 ***
## DHW4                 0.12225    0.04909   2.490 0.023399 *  
## recent.meanDHW      -0.16296    0.08962  -1.818 0.086680 .  
## returnDHW6           0.05856    0.02546   2.300 0.034371 *  
## SpeciesP. meandrina  0.41367    0.14079   2.938 0.009184 ** 
## SpeciesP. verrucosa  0.68668    0.16006   4.290 0.000495 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2882 on 17 degrees of freedom
## Multiple R-squared:  0.7835, Adjusted R-squared:  0.7199 
## F-statistic: 12.31 on 5 and 17 DF,  p-value: 3.718e-05
ggeffects::ggeffect(lm.dredge1) %>% plot
## $DHW4

## 
## $recent.meanDHW

## 
## $returnDHW6

## 
## $Species

The best model is ED50 ~ DHW4 + recent.meanDHW + returnDHW6 + Species Only varSST was not included. However, we need to also consider if the species interaction is significant.

Dredging with species interaction.

Dredging across all selected variables suggest a model that includes DHW4, MMM and recent.meanDHW best represents the data so far. However, we’ve seen there may be an interaction between Species * DHW4. Before forcing this interaction using dredge, we will see which model yields best AIC score with interaction and additional fixed effects.

lm.1 = lm(ED50~Species + DHW4 + recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.2 = lm(ED50~Species * DHW4 + recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.3 = lm(ED50~Species * returnDHW6 + DHW4 + recent.meanDHW, data = dredge.dat1)
lm.4 = lm(ED50~Species * DHW4 * returnDHW6 + recent.meanDHW, data = dredge.dat1)


AICc(lm.1, lm.2, lm.3, lm.4)
##      df     AICc
## lm.1  7 22.56228
## lm.2  9 21.94864
## lm.3  9 28.24851
## lm.4 14 61.85056
#lm.2 has the lowest AIC, suggesting it is the better model. 
anova(lm.1, lm.2)
## Analysis of Variance Table
## 
## Model 1: ED50 ~ Species + DHW4 + recent.meanDHW + returnDHW6
## Model 2: ED50 ~ Species * DHW4 + recent.meanDHW + returnDHW6
##   Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
## 1     17 1.41235                              
## 2     15 0.87572  2   0.53663 4.5959 0.02774 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#The two models are also sufficiently different from each other to include the interaction
summary(lm.2)
## 
## Call:
## lm(formula = ED50 ~ Species * DHW4 + recent.meanDHW + returnDHW6, 
##     data = dredge.dat1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38367 -0.08972  0.03304  0.09462  0.43907 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.78952    0.76071   7.611 1.58e-06 ***
## SpeciesP. meandrina       1.84012    0.55064   3.342 0.004460 ** 
## SpeciesP. verrucosa       2.19242    0.60783   3.607 0.002589 ** 
## DHW4                      0.25490    0.06028   4.229 0.000729 ***
## recent.meanDHW           -0.17568    0.07644  -2.298 0.036344 *  
## returnDHW6                0.05255    0.02162   2.431 0.028076 *  
## SpeciesP. meandrina:DHW4 -0.20679    0.07812  -2.647 0.018295 *  
## SpeciesP. verrucosa:DHW4 -0.21270    0.08280  -2.569 0.021379 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2416 on 15 degrees of freedom
## Multiple R-squared:  0.8658, Adjusted R-squared:  0.8031 
## F-statistic: 13.82 on 7 and 15 DF,  p-value: 1.633e-05
summary(lm.3)
## 
## Call:
## lm(formula = ED50 ~ Species * returnDHW6 + DHW4 + recent.meanDHW, 
##     data = dredge.dat1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43125 -0.12992  0.00647  0.14245  0.57188 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     6.71358    0.80810   8.308 5.39e-07 ***
## SpeciesP. meandrina            -0.08909    0.30582  -0.291   0.7748    
## SpeciesP. verrucosa             0.35133    0.37683   0.932   0.3659    
## returnDHW6                      0.01414    0.03497   0.404   0.6917    
## DHW4                            0.13149    0.04793   2.744   0.0151 *  
## recent.meanDHW                 -0.14462    0.08738  -1.655   0.1187    
## SpeciesP. meandrina:returnDHW6  0.08080    0.04397   1.838   0.0860 .  
## SpeciesP. verrucosa:returnDHW6  0.05375    0.05215   1.031   0.3190    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2771 on 15 degrees of freedom
## Multiple R-squared:  0.8235, Adjusted R-squared:  0.7411 
## F-statistic: 9.997 on 7 and 15 DF,  p-value: 0.0001144
#the interaction with DHW4 is important but not the interaction with returnDHW6 so will be used in the full dredge model.

The interaction is significant and the two models are significantly different from one another so we should include the interaction. (We also explored a three-way interaction with MMM and the interaction between Species and MMM, neither were significant)

vif(lm.2)
##                      GVIF Df GVIF^(1/(2*Df))
## Species        482.320105  2        4.686341
## DHW4             3.455627  1        1.858932
## recent.meanDHW   2.111638  1        1.453148
## returnDHW6       1.734821  1        1.317126
## Species:DHW4   517.710799  2        4.770038

In addition there is no colinearity between these core variables, which is good.

dredge.mod2 <- dredge(lm(ED50~Species*DHW4 +., data = dredge.dat1), rank = "AIC", m.lim = c(1,4))
## Fixed term is "(Intercept)"
#just look at the top scoring models 
head(dredge.mod2, 10)
## Global model call: lm(formula = ED50 ~ Species * DHW4 + ., data = dredge.dat1)
## ---
## Model selection table 
##    (Int)     DHW rcn.mDH     rDH Spc    vSS DHW:Spc df logLik  AIC delta weight
## 46 4.321 0.30680         0.08390   +              +  8  1.479 13.0  0.00  0.391
## 44 6.986 0.21650 -0.2929           +              +  8  1.129 13.7  0.70  0.275
## 16 6.607 0.12220 -0.1630 0.05856   +                 7 -0.548 15.1  2.05  0.140
## 58 1.724 0.22950                   + 1.4590       +  8 -0.365 16.7  3.69  0.062
## 14 5.212 0.17540         0.08732   +                 6 -2.592 17.2  4.14  0.049
## 30 4.335 0.14940         0.06693   + 0.4729          7 -2.277 18.6  5.51  0.025
## 12 7.998 0.06906 -0.2910           +                 6 -3.664 19.3  6.29  0.017
## 28 5.601 0.06648 -0.1862           + 0.7482          7 -2.720 19.4  6.40  0.016
## 11 8.734         -0.3401           +                 5 -4.907 19.8  6.77  0.013
## 27 6.183         -0.2278           + 0.7871          6 -3.967 19.9  6.89  0.012
## Models ranked by AIC(x)
#extracting the best fit model using AIC score 
bestmodel <- get.models(dredge.mod2, 1)[[1]]
lm.dredge2 <- lm(bestmodel, data = dredge.dat1)
summary(lm.dredge2)
## 
## Call:
## lm(formula = bestmodel, data = dredge.dat1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.46658 -0.17588  0.03949  0.13706  0.58307 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               4.32105    0.46482   9.296 7.51e-08 ***
## DHW4                      0.30680    0.06293   4.876 0.000168 ***
## returnDHW6                0.08390    0.01889   4.442 0.000410 ***
## SpeciesP. meandrina       1.90207    0.61922   3.072 0.007300 ** 
## SpeciesP. verrucosa       1.94108    0.67319   2.883 0.010806 *  
## DHW4:SpeciesP. meandrina -0.21359    0.08789  -2.430 0.027230 *  
## DHW4:SpeciesP. verrucosa -0.18608    0.09230  -2.016 0.060924 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.272 on 16 degrees of freedom
## Multiple R-squared:  0.8185, Adjusted R-squared:  0.7504 
## F-statistic: 12.03 on 6 and 16 DF,  p-value: 3.728e-05
ggeffects::ggeffect(lm.dredge2) %>% plot
## $DHW4

## 
## $returnDHW6

## 
## $Species

When we now include the interaction recent.meanDHW drops out from the best model. However, it is still important and results in a better model.

lm.2 = lm(ED50~Species * DHW4 + recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.5 = lm(ED50~Species * DHW4 + returnDHW6, data = dredge.dat1)
AICc(lm.2, lm.5)
##      df     AICc
## lm.2  9 21.94864
## lm.5  8 23.32708
anova(lm.2, lm.5)
## Analysis of Variance Table
## 
## Model 1: ED50 ~ Species * DHW4 + recent.meanDHW + returnDHW6
## Model 2: ED50 ~ Species * DHW4 + returnDHW6
##   Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
## 1     15 0.87572                              
## 2     16 1.18410 -1  -0.30838 5.2821 0.03634 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple interactions

lm.2 = lm(ED50~Species * DHW4 + recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.6 = lm(ED50~Species * DHW4 * recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.7 = lm(ED50~Species * DHW4 * returnDHW6 + recent.meanDHW, data = dredge.dat1)
lm.8 = lm(ED50~Species * returnDHW6 + DHW4  + recent.meanDHW, data = dredge.dat1)
lm.9 = lm(ED50~Species * recent.meanDHW + DHW4  + returnDHW6, data = dredge.dat1)

AICc(lm.2, lm.6, lm.7, lm.8, lm.9)
##      df     AICc
## lm.2  9 21.94864
## lm.6 14 66.28199
## lm.7 14 61.85056
## lm.8  9 28.24851
## lm.9  9 29.57824

None of the other interactions improved the model.

Model validation

library(performance)
check_model(lm.2)
## Loading required namespace: qqplotr
## For confidence bands, please install `qqplotr`.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 22 rows containing missing values (geom_text_repel).

r.squaredGLMM(lm.2)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.8147412 0.8147412
summary(lm.2)
## 
## Call:
## lm(formula = ED50 ~ Species * DHW4 + recent.meanDHW + returnDHW6, 
##     data = dredge.dat1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38367 -0.08972  0.03304  0.09462  0.43907 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.78952    0.76071   7.611 1.58e-06 ***
## SpeciesP. meandrina       1.84012    0.55064   3.342 0.004460 ** 
## SpeciesP. verrucosa       2.19242    0.60783   3.607 0.002589 ** 
## DHW4                      0.25490    0.06028   4.229 0.000729 ***
## recent.meanDHW           -0.17568    0.07644  -2.298 0.036344 *  
## returnDHW6                0.05255    0.02162   2.431 0.028076 *  
## SpeciesP. meandrina:DHW4 -0.20679    0.07812  -2.647 0.018295 *  
## SpeciesP. verrucosa:DHW4 -0.21270    0.08280  -2.569 0.021379 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2416 on 15 degrees of freedom
## Multiple R-squared:  0.8658, Adjusted R-squared:  0.8031 
## F-statistic: 13.82 on 7 and 15 DF,  p-value: 1.633e-05
#Slopes comparison: 
emtrends(lm.2, pairwise~Species, var = "DHW4")
## $emtrends
##  Species       DHW4.trend     SE df lower.CL upper.CL
##  A. cf humilis     0.2549 0.0603 15   0.1264    0.383
##  P. meandrina      0.0481 0.0581 15  -0.0758    0.172
##  P. verrucosa      0.0422 0.0693 15  -0.1054    0.190
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                     estimate     SE df t.ratio p.value
##  A. cf humilis - P. meandrina  0.20679 0.0781 15 2.647   0.0455 
##  A. cf humilis - P. verrucosa  0.21270 0.0828 15 2.569   0.0528 
##  P. meandrina - P. verrucosa   0.00591 0.0821 15 0.072   0.9971 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
library(jtools)
library(interactions)
## 
## Attaching package: 'interactions'
## The following objects are masked from 'package:jtools':
## 
##     cat_plot, interact_plot, johnson_neyman, probe_interaction,
##     sim_slopes
sim_slopes(lm.2, pred = DHW4, modx = Species, johnson_neyman = FALSE)
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of DHW4 when Species = P. verrucosa: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   0.04   0.07     0.61   0.55
## 
## Slope of DHW4 when Species = P. meandrina: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   0.05   0.06     0.83   0.42
## 
## Slope of DHW4 when Species = A. cf humilis: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   0.25   0.06     4.23   0.00
emtrends(lm.2, ~Species, var = "recent.meanDHW")
##  Species       recent.meanDHW.trend     SE df lower.CL upper.CL
##  A. cf humilis               -0.176 0.0764 15   -0.339  -0.0128
##  P. meandrina                -0.176 0.0764 15   -0.339  -0.0128
##  P. verrucosa                -0.176 0.0764 15   -0.339  -0.0128
## 
## Confidence level used: 0.95
sim_slopes(lm.2, pred = recent.meanDHW, modx = Species, johnson_neyman = FALSE)
## Warning: recent.meanDHW and Species are not included in an interaction with one
## another in the model.
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of recent.meanDHW when Species = P. verrucosa: 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.18   0.08    -2.30   0.04
## 
## Slope of recent.meanDHW when Species = P. meandrina: 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.18   0.08    -2.30   0.04
## 
## Slope of recent.meanDHW when Species = A. cf humilis: 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.18   0.08    -2.30   0.04

###Relative Importance of Predictors

library('relaimpo')
#calc.relimp(lm.2, type = "lmg")
calc.relimp(lm.2, type = "lmg")@lmg 
##        Species   Species:DHW4           DHW4 recent.meanDHW     returnDHW6 
##     0.28829712     0.08514947     0.16291827     0.18219391     0.14721637
sum(calc.relimp(lm.2, type = "lmg")@lmg)
## [1] 0.8657751
calc.relimp(lm.2, type = "lmg")@lmg  /  sum(calc.relimp(lm.2, type = "lmg")@lmg)
##        Species   Species:DHW4           DHW4 recent.meanDHW     returnDHW6 
##     0.33299307     0.09835056     0.18817619     0.21044022     0.17003996

For the output, we see that total proportion of variance explained by model is 86.6%. Species contributed 28.8% of variance, interaction of species * DHW4 contributed 8.5%, DHW4 contributed 16.3%, recent mean DHW contributed 18.2% and returnDHW6 accounted for 14.7% (total of these predictors = 86.6%, which is the total variance explained by the model)

Model Predictions

Fig3a =  ggplot(mod.sep, aes(y = emmean, x = Reef, ymin = lower.CL, ymax = upper.CL, col = Species, shape = Species)) +
  geom_pointrange(position = position_dodge(width = .4), size = .2) +
  scale_y_continuous(expand = c(0,0), lim = c(6.1,8.6), breaks = c(6.5, 7, 7.5, 8, 8.5)) +
  scale_colour_manual(values=col) +
  labs(x="", 
       y="ED50 (95% CI)") +
  theme_classic() + 
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        legend.title = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        legend.position = "bottom",
        legend.direction = "vertical",
        legend.text = element_text(face = "italic", size = 8, family = "Helvetica"),
        legend.key.size = unit(3, units = "mm"),
        legend.margin = margin(0,0,0,0, unit = "mm"),
        axis.text = element_text(size = 8, family = "Helvetica"),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title = element_text(size = 9, family = "Helvetica"),)

Fig3a

ggsave("2_figures/Fig3_ED50.pdf", plot = Fig3a, width = 70, height = 78, units = "mm")
#DHW4
grid.DHW4 = with(dredge.dat1, list(DHW4 = seq(min(DHW4)-.1, max(DHW4)+.1, len = 50)))
em.DHW4 = emmeans(lm.2, ~DHW4|Species, at = grid.DHW4, type = "response") %>% 
  as.data.frame() %>% 
  mutate(var = "DHW4") %>% 
  rename(value = DHW4)
#returnDHW6
grid.returnDHW6 = with(dredge.dat1, list(returnDHW6 = seq(min(returnDHW6)-.1, max(returnDHW6)+.1, len = 50)))
em.returnDHW6 = emmeans(lm.2, ~returnDHW6|Species, at = grid.returnDHW6, type = "response") %>% 
  as.data.frame() %>% 
  mutate(var = "returnDHW6") %>% 
  rename(value = returnDHW6)
#recent.meanDHW
grid.recent.meanDHW = with(dredge.dat1, list(recent.meanDHW = seq(min(recent.meanDHW)-.1,
                                                                  max(recent.meanDHW)+.1, len = 50)))
em.recent.meanDHW = emmeans(lm.2, ~recent.meanDHW|Species, at = grid.recent.meanDHW, type = "response") %>% 
  as.data.frame() %>% 
  mutate(var = "recent.meanDHW") %>% 
  rename(value = recent.meanDHW)

library(wesanderson)
col = wes_palette("GrandBudapest1", 4, type = "continuous")

DHW4.lm.plot = ggplot(em.DHW4, aes(x = value, y = emmean, fill = Species)) +
  geom_point(data = dredge.dat1, aes(y = ED50, x = DHW4, col = Species, shape = Species), alpha = .8) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = .2) +
  geom_line(aes(col = Species)) +
  scale_color_manual(values = col) +
  scale_fill_manual(values = col) +
  scale_y_continuous(expand = c(0,0), lim = c(6.2,8.6), breaks = c(6.5, 7, 7.5, 8, 8.5)) +
  scale_x_continuous(expand = c(0,0), lim = c(5,11), breaks = c(6,  8, 10)) +
  coord_cartesian(xlim = c(5.6,10.6), ylim = c(6.1,8.6)) +
  labs(y = "ED-50 (95% CI)", x = "# events > 4°C-weeks") +
  theme_classic() +
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        legend.position = "none", 
        axis.text = element_text(size = 8, family = "Helvetica"),
        axis.title = element_text(size = 9, family = "Helvetica"),
        axis.title.y = element_blank()) 

returnDHW6.lm.plot = ggplot(em.returnDHW6, aes(x = value, y = emmean, fill = Species)) +
  geom_point(data = dredge.dat1, aes(y = ED50, x = returnDHW6, col = Species, shape = Species), alpha = .8) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = .2) +
  geom_line(aes(col = Species)) +
  scale_color_manual(values = col) +
  scale_fill_manual(values = col) +
  scale_y_continuous(expand = c(0,0), breaks = c(6.5, 7, 7.5, 8, 8.5)) + #lim = c(6.2,8.6)
  scale_x_continuous(expand = c(0,0), breaks = c(3,6, 9, 12)) + 
  coord_cartesian(xlim = c(1.85, 13.09), ylim = c(6.1,8.6)) +
  labs(y = "ED-50 (95% CI)", x = "Return time events > 6°C-weeks") +
  theme_classic() +
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        legend.position = "none", 
        axis.text = element_text(size = 8, family = "Helvetica"),
        axis.title = element_text(size = 9, family = "Helvetica"),
        axis.title.y = element_blank())

recent.meanDHW.lm.plot = ggplot(em.recent.meanDHW, aes(x = value, y = emmean, fill = Species)) +
  geom_point(data = dredge.dat1, aes(y = ED50, x = recent.meanDHW, col = Species, shape = Species), alpha = .8) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = .2) +
  geom_line(aes(col = Species)) +
  scale_color_manual(values = col) +
  scale_fill_manual(values = col) +
  scale_y_continuous(expand = c(0,0), breaks = c(6.5, 7, 7.5, 8, 8.5)) + #lim = c(6.2,8.6)
  scale_x_continuous(expand = c(0,0), breaks = c(4,5,6)) + #lim = c(27.2,28.8), 
  coord_cartesian(xlim = c(3.44,6.6), ylim = c(6.1,8.6)) +
  labs(y = "ED-50 (95% CI)", x = "mean maxDHW 2016-20") +
  theme_classic() +
  theme(axis.line = element_blank(),
        panel.border = element_rect(size = .5, fill = "transparent"),
        legend.position = "none", 
        axis.text = element_text(size = 8, family = "Helvetica"),
        axis.title = element_text(size = 9, family = "Helvetica"),
        axis.title.y = element_blank()) 
gpt = theme(axis.title = element_blank(), 
      plot.margin = unit(c(1,1,1,1), "mm"))

ggarrange(Fig3a + theme(legend.position = "none"),
          DHW4.lm.plot + gpt,
          returnDHW6.lm.plot + gpt,
          recent.meanDHW.lm.plot + gpt,
          ncol = 4, align = "hv")

ggsave("2_figures/Fig3_model_v3.pdf", width = 180, height = 53, units = "mm")